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Regarding  ongoing  roaaarch,  noat  efforts  have  been  devoted  to  2D  FD 
synthetics  for  siaulating  seismic  save  propagation  in  a  complex  lithosphere.  In 
practice  model  set-up  for  FO  computations  have  proved  cumbersome  so  se  are 
elaborating  on  automating  structural  model  generation.  The  nes  approach  is  tied 
to  using  an  optical  scanner  for  reading  a  figure  of  the  multilayered  crustal 
model. 

Event  magnitude  is  an  versatile  parameter  in  many  seismological  studies  not 
at  least  in  the  context  of  nuclear  test  ban  monitoring.  For  simulating  netvork 
signal  detection  capabilities  the  noise  level  at  individual  stations  are  needed 
srhich  in  practice  appears  difficult  to  obtain.  Ideally,  we  uould  prefer  to  have 
such  noise  level  estimates  available  on-line  to  measure  the  ’on-line  threshold 
monitoring’  performance.  Most  signal  detections  in  use  are  of  the  STA/LTA  type 
where  the  short-term-average  (STA)  is  in  fact  an  a  RMS  trace  estimate.  Initial 
test  attempting  to  relate  STAs  to  maximum  amplitude  used  in  local  magnitude 
formulas  shows  that  error  less  than  0.1  magnitude  unit  should  be  feasible. 

The  new  Sonseca  (ESLA)  seismic  array  near  Toledo.  Spain,  has  recently  become 
operational.  For  its  operation,  crustal  structure  models  are  of  some 
importance.  Using  Pn  and  Sn  travel  time  observations  from  the  Iberian 
seismograph  network  we  have  derived  a  2D  map  of  subMoho  P-  and  S-velocity 
variations  for  this  area.  The  central  parts  of  the  region  appears  to  be  rather 
homogeneous  while  the  largest  anomalies  are  to  the  south.  Crustal  thickness 
variations  are  pronounced  ranging  from  about  20  to  40  km. 

Finally,  a  comprehensive  study  aimed  at  seismic  event  classification  at 
local  and  regional  distances  has  been  starded.  Preparatory  efforts  here 
comprise  software  development  and  selection  of  a  comprehensive  data  base.  The 
research  strategy  is  aimed  at  examining  the  discriminative  potential  of  any 
part  of  the  seismic  records  as  bracketed  by  the  group  velocity  interval  8.5  to 
2.5  km/sec.  Special  attention  would  be  given  to  focal  depth  estimation  through 
inversion  of  relatively  low-frequency  recordings  for  small  events  at  local  and 
regional  distances. 


Unclassified 


Contents 


1  Ovwrvww  2 

1.1  Objectives .  2 

1.2  Summary  of  ongoing  research .  2 

1.3  Summaries  of  studies  completed .  3 

1 .4  Outline  of  the  report .  3 

2  Ongoing  research  activities  5 

2.1  Seismic  event  classification .  3 

2.2  Real-Time  Magnitude  Elstimation .  8 

2.3  Crustal  structure  of  Iberia .  12 

2.4  Crustal  modelling  for  2D  synthetics .  20 

3  Studies  completed  24 

3.1  Spatial  and  temporal  patterns  of  the  Fennoscandian  seismicity  -  an  exercise  in 

explosion  monitoring .  24 

3.2  2D  finite  difference  elastic  wave  modeling  including  surface  topography .  25 

3.3  Seismic  wave  propagation  in  complex  crust  -  upper  mantle  media  using  2D  finite 

difference  synthetics .  26 

4  Future  work  plans  27 

5  Appendix  A  A>1 

6  Appendix  B  B-1 


•  • 


•  • 


•  • 


•  •  • 


•  • 


DTIC  QUAUTY  araPBCTED  s 


Aecesion  t^cr  ^  — 

NTIS 

CRA4I  U 

DTIC 

T.A3  r; 

U'ia.'i 

•iji:  ;•  i'-r' 

Oi-f  1 

.•t.o-/ 

1  «va...ioi  =fy  i 

Dist 

i  A  VO.; 

/l-l 

i 

•  •  • 


1  Overview 

1.1  Objectives 

•  Foster  a  better  understanding  of  seismic  wave  propagation  in  complex  lithospheric  media. 

•  Principal  modelling  tool  is  2D  finite  difference  (FD)  synthetics  for  models  also  including 
surface  topography. 

•  Corroborating  synthetics  through  analysis  of  array  and  network  data. 

•  Use  ’synthetic’  knowledge  in  design  of  event  classification  schemes. 

•  Seismic  network  performance. 

Note,  a  basic  element  in  research  is  curiosity  which  in  practice  will  imply  a  certain  flexibility 
regaiding  adherence  to  objectives. 

1.2  Summary  of  ongoing  research 

Most  efforts  have  been  devotes  to  2D  FD  synthetics  for  simulating  seismic  wave  propagation 
in  a  complex  lithosphere.  The  FD  scheme  used  is  tied  to  the  numerical  solution  of  the  elastody- 
namic  wave  equation.  Spatial  partial  differentiation  is  achieved  through  cost-optimized,  dispersion- 
bounded,  high-order  finite  difference  operators  on  a  staggered  grid.  For  time  stepping,  a  leap-frog 
technique  is  used.  For  the  numerical  dispersion  relations,  the  stability  limit  and  btindwidth  intro¬ 
duced  by  the  discretization,  the  reader  is  referred  to  Sguazzero,  Kindeland  and  Kamel  (1989). 

Recently  we  have  demonstrated  the  usefulness  of  this  approach  by  computing  synthetic  seismo¬ 
grams  for  a  set  of  progressively  complex  lithosphere  models  (Hestholm,  Husebye  and  Ruud.  1994). 
In  practice,  model  set-up  for  FD  computations  have  proved  cumbersome  so  we  are  elaborating  on 
automating  structural  model  generation.  The  new  approach  is  tied  to  using  a  scttnner  for  reading  a 
published  multilayered  crustal  model.  For  any  model  point  proper  layer  assignment  is  performed 
automatically  and  hence  a  corresponding  specification  of  physical  parameters  like  velocity  and 
density.  The  computed  synthetics  are  SEG-Y  formatted  and  thus  can  easily  be  processed  using 
standard  softweue  packages  like  the  ProMAX. 

Event  magnitude  is  tm  versatile  parameter  in  many  seismologicai  studies  not  at  least  in  the 
context  of  nuclear  test  ban  monitoring.  For  simulating  network  signal  detection  capabilities  the 
noise  level  at  individual  stations  are  needed  which  in  practice  appears  to  be  difficult  to  obtain. 
Ideally,  we  would  prefer  to  have  such  noise  level  estimates  available  on-line  to  measure  the  ’on-line 
threshold  monitoring’  performance  (Ringdal  and  Kvaema,  1992).  To  achieve  this  we  are  exploring 
the  usefulness  of  a  random  vibration  theory  formula,  namely  (Cartwright  and  Longnet- Higgins, 
1956). 

^mox  “  f{N)A  rm$  (1) 

where  Amax  is  max  amplitude;  f(N)  where  N  is  the  number  of  extrema  with  the  RMS-window 
tmd  Artrn  is  amplitude  root-mean-square.  Most  signal  detections  in  use  are  of  the  STA/LTA  type 
where  the  short-term-average  (STA)  is  in  fact  an  a  (RMS)  trace  estimate.  Initial  test  results  from 
real  recordings  indicate  that  the  above  formula  is  valid;  error  estimates  amount  to  ±0.1  magnitude 
unit.  In  other  word,  on-line  magnitude  estimation  should  be  feasible. 


Tlie  aew  Sonseca  (ESLA)  seismic  array  near  Toledo,  Spain,  has  become  operational.  For  its 
operation,  crustal  structure  models  are  of  some  importance  for  slowness-vector  calibration  and 
related  purposes.  Using  Pn  auid  Sn  travel  time  observations  from  the  Iberian  seismograph  network 
with  supplementary  observations  &om  stations  in  France,  Algeria  and  Morocco  we  have  derived  a 
2D  map  of  subMoho  P-  and  S-velocity  variations  for  Iberia  and  adjacent  areas.  The  central  parts 
of  the  region  appears  to  be  rather  homogeneous  while  the  largest  anomalies  are  found  in  the  Betics 
and  the  Alboran  Sea  to  the  south.  Crustal  thickness  variations  are  also  pronounced  in  this  area 
ranging  from  c.  20  to  40  km. 

Finally,  a  comprehensive  study  aimed  at  seismic  event  classification  at  local  and  regional  dis¬ 
tances.  Preparatory  efforts  here  comprise  software  development  and  selection  of  a  comprehensive 
data  base;  we  opted  for  the  NORESS  and  ARCESS  events  used  by  Lacoss  et  al.  (1992).  The 
research  strategy  is  turned  at  examining  the  discriminative  potential  of  any  part  of  the  seismic 
records  as  bracketed  by  the  group  velocity  interval  8.5  to  2.5  km/sec.  Special  attention  would  be 
given  to  focal  depth  estimation  through  inversion  of  relatively  low-  frequency  recordings  for  small 
events  at  local  and  regional  distances.  Various  schemes  for  suppression  of  propagating  noise  on  3 
component  recordings  are  now  under  considerations;  if  successful  we  experiment  with  wave  form 
inversion  for  accurate  estimation  of  focal  depths. 

1.3  Summaries  of  studies  completed 

In  addition  to  the  paper  (Hestholm  et  al,  1994)  included  in  the  First  Scientific  Report  submitted 
to  AFOSR  (  19.  May  1993)  two  other  studies  have  been  completed  in  the  reporting  period:  Tar- 
vainen  and  Husebye  (1994)  and  Hestholm  and  Ruud  (1994).  The  first  of  these  reports  deab  with 
spatial  and  temporal  patterns  in  Fennoscandian  seismicity  essentially  aunounting  to  an  exercise 
in  monitoring  mining  amd  other  explosions.  We  demonstrate  that  there  exist  a  very  consistent 
diurnal/weekly  and  spatial  patterns  for  the  numerous  e.xplo6ion  events  included  in  local  seismolog- 
ical  bulletins.  The  seismological  epicenter  locations  are  clearly  inferior  to  these  made  by  a  skilled 
analyst  and  tied  to  waveform  stationarity  for  a  given  mine. 

The  other  study  (Hestholm  and  Ruud,  1994)  presents  a  numerical  2D  finite  difference  (FD) 
algorithm  for  including  free  surface  topography  in  synthetic  seismogram  2tnalysis.  The  realism  of 
this  report  is  also  demonstrated  through  displays  of  snapshots  of  near  surface  effects  for  models 
without  and  with  topography  included. 

1.4  Outline  of  the  report 

The  First  Technical  Report  under  Contract  AFOSR  Grant  F49620-92-5-0510  is  divided  in  2  parts; 
in  Section  2  we  give  details  on  ongoing  research  while  in  Section  3  we  present  abstracts  of  studies 
completed  suid  also  discuss  some  implications  of  the  results  obtained  here.  The  completed  text  of 
these  2  reports  are  given  in  Appendix  A  and  B. 
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2  Ongoing  research  activities 


Here  w«  give  details  on  ongoing  research  activities  some  of  which  are  scheduled  for  completion 
at  15.  February  1994  i.e.,  the  deadline  for  contributions  for  the  proceedings  of  the  WORKSHOP 
ON  PLANNING  AND  PROCEDURES  FORGSETT-3  held  in  Erice,  Italy  11-13  November  1993 
As  is  well  known  GSEnT-3  is  short  for  group  of  Scientific  Expert  Technical  Test  no.  3  which 
is  an  exercise  in  testing  various  aspects  of  a  prototype  global  seismic  monitoring  system  for  a 
potential  comprehensive  nuclear  test  treaty  (CTBT).  The  GSE  reports  to  the  UN  Conference  on 
Disarmament  in  Geneva. 

2.1  Seismic  event  classification 

In  recent  years  there  has  been  renewed  interest  in  event  classification  of  small  local  events  which 
remain  problematic.  The  many  studies  undertaken  here  reflect  the  renewed  political  interest  in  a 
comprehensive  test  ban  treaty  effectively  prohibiting  any  test  of  nuclear  weapons.  Most  promising 
discriminant  at  locsd  disttmces  appears  to  be  the  spectral  ratio  of  Pn/Lg-phases  in  the  frequency 
band  2-3  Hz.  Physically  the  explanation  is  that  earthquakes  source  excitation  is  relative  effective 
in  generating  shear  waves.  Since  up  to  now  most  nuclear  testing  are  confined  to  a  few  sites,  a 
fair  amount  of  classification  research  is  tied  to  discriminating  between  earthquakes  and  chemical 
explosions  stemming  from  mining  and  quarry  blasting  activities.  In  comparison  to  nuclear  tests 
the  chemical  explosions  are  truly  numerous  so  the  challenge  is  naturally  to  be  able  to  differentiate 
between  the  3  major  types  of  seismic  sources,  namely  earthquakes  (EQ),  nuclear  explosions  (NE) 
and  chemical  explosions  (CE).  For  convenience  the  explosions  are  referred  to  as  EX  unless  we 
need  to  distinguish  between  NE  and  CE.  A  distinct  signal  attribute  of  CE  is  the  so  called  spectral 
scalloping  effect  caused  by  the  commonly  used  technique  of  ripple  firing  typicsd  of  lauge  chemical 
explosions  (Hedlin  et  al,  1990). 

The  recent  trend  in  seismic  event  classification  research  is  to  introduce  a  multitude  of  potential 
discriminants  and  then  use  non-statistical  techniques  like  artificial  neural  networks  (ANN)  schemes 
to  ensure  optimum  performances.  A  possible  drawback  with  the  ANN  approach  is  that  the  learning 
set  has  to  be  large  (say  30  events  each  for  the  EQ  and  EX  populations)  and  also  that  there  is 
no  easy  assessment  of  the  relative  merits  of  the  individual  discriminants.  In  contrast,  statisticed 
approaches  like  those  based  on  pattern  recognition  schemes  (Tjostheim,  1982;  Tsvang  et  al.  1993) 
function  well  for  relative  small  event  populations  and  also  enable  us  to  estimate  the  relative 
significance  of  individual  discriminants.  The  latter  technique  can  be  extended  to  handle  one-sided 
event  populations  that  is  explosions  only  or  earthquakes  only.  For  strong  events,  recorded  by 
many  stations  in  the  teleseismic  window,  seismic  event  classification  is  not  problematic  and  even 
traditional  mb.Ms  and  complexity  discriminants  are  very  efficient. 

Seismic  Discriminant  Features 

Numerous  studies  on  seismic  event  classification  during  the  last  decade  have  been  based  on 
spectrtd  characteristics  of  crust  and  and  upper  mantle  phases  like  Pn,  Pg,  Lg,  Rg  etc.  or  on  their 
spectral  ratios  in  a  multitude  of  frequency  bands.  This  in  turn  reflect  extensive  observational 
experiences  that  explosions  and  earthquakes  produce  different  spectral  ^unplitudes  for  the  same 
direct  phase  at  comparable  distance  ranges  and  wave  paths  (Pulli  and  Dysart,  1992).  A  potential 
drawback  with  this  class  of  discriminrmts  is  that  the  direct  wave  phases  can  be  strongly  affected 


% 


m 


•  • 


•  • 


•  •  • 


•  • 


5 


by  a  particuiar  source-receiver  path  thus  making  it  difficult  to  separate  source  and  propagation 
effects.  In  the  extreme  a  certain  discriminant  may  fail  when  the  EXs  and  the  EQs  are  not  closely 
located.  A  related  problem  is  that  late  P-coda  waves,  direct  Sn,  Lg  and  surface  waves  are  mixed 
together  and  the  isolation  of  separate  phases  are  not  trivial.  In  this  context,  Lg-coda  waves  exhibit 
a  distinct  advantage  as  a  characteristic  feature  here  is  that  their  amplitude  decay  with  distance 
reflects  the  average  crustal  properties  of  the  region  encompassing  the  source  and  the  receivers 
(e.g.,  see  Su  et  al.  1991).  In  other  words  we  can  separate  source  and  path  effects  on  the  coda 
waves.  Likewise,  the  spectral  scalloping  effect  stemming  from  ripple  firing  in  mining  operations 
appears  to  be  most  distinctive  in  the  coda  (Hedlin  et  al.,  1990). 

Data  Analysts  Strategy 

The  data  at  hand  for  classification  analysis  stem  from  array  and  3-component  (3C)  station 
recordings  from  events  at  local  and  regional  distance  ranges.  This  in  turn  would  entml  that  a 
flexible  analysis  strategy  should  be  implemented  simply  because  signal  and  coda  durations  are 
distance  dependent.  Another  problem  is  that  of  introducing  corrections  for  geometrical  spreading 
and  attenuation  both  of  which  appear  to  vary  considerably  from  one  region  to  another  and  also 
within  a  given  region.  Such  corrections,  often  ignored  in  many  event  classification  studies,  should 
be  introduced  as  in  case  of  the  powerful  Pn/Lg  discriminant  both  the  geometricad  spreading  and 
attenuation  are  different  for  these  two  phases.  Linear  event  scaling  would  also  be  important  except 
for  spectral  ratio  type  of  discriminants  because  we  undertake  a  comparison  of  event  with  different 
magnitudes.  With  array  data  at  hand  we  want  to  incorporate  measures  to  facilitate  assessment 
of  the  relative  merits  of  such  recording  systems  vis-a-vis  modern  3C  stations.  Finally,  to  rettiin 
flexibility  we  decided  from  the  very  onset  to  perform  our  analysis  on  sonogram  data  (signal  power 
as  a  function  of  time  and  frequency  averaged  over  the  array  elements)  with  the  added  advantage 
of  reducing  the  storage  requirements  of  the  original  time  domain  recordings. 

We  recognize  the  following  principal  part  of  any  local/regional  seismic  record;  Preceeding  noise, 
Pn,  Pg,  P-coda,  Sn,  Lg  and  Lg-coda.  The  noise  is  of  duration  30  to  60  sec  depending  of  available 
record  length.  The  body  wave  phases  are  of  length  3.2  to  6.4  sec  and  distance  dependent.  The 
P-coda  is  the  segment  between  Pn/Pg-end  and  Sn/Sg-start.  The  Lg  is  taken  to  start  a  sec  prior 
to  given  start  time  or  alternatively  to  the  equivalent  group  velocity  onset  time  of  3.6  km/sec. 
Lg  ends  at  the  ’group  velocity’  time  of  3.2  km/sec.  The  Lg-coda  ends  at  an  equivalent  group 
velocity  of  2.0  km/sec.  We  will  experiment  with  a  coda  subdivision  in  3  segments  that  is  of  group 
velocity  lengths  of  3. 2-2. 8;  2.8-2.4  and  2. 4-2.0  km/sec  respectively.  The  basic  spectrum  analysis 
is  performed  via  the  Fast  Fourier  Transform  (FFT)  for  data  windows  of  unit  length  3.2  sec  (128 
samples)  with  Hanning  type  tapering.  The  spectrum  is  updated  at  1.6  sec  intervals,  and  additional 
smoothing  is  introduced  by  stacking  of  spectrum  segments. 

Study  Summary 

So  far  most  of  our  efforts  have  been  tied  to  software  development  and  getting  access  to  a  large 
data  base  of  seismic  array  events.  It  is  not  quite  trivial  to  use  these  excellent  CSS-generated  data 
bases  but  practical  advises  &om  R.  Lacoss  (Lincoln  Lab.)  proved  very  helpful.  Adapting  to  a  new 
computer  environment  in  Bergen  have  also  taken  time.  Anyway,  we  are  now  in  the  position  to 
start  large  scale  discriminant  feature  extraction  from  a  large  number  of  local  and  regional  events 
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as  recorded  by  the  NORESS  and  ARCESS  arrays. 
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2.2  Real-Time  Magnitude  Estimation 

The  problem  of  estimating  the  size  of  an  earthquake  (EQ)  has  intrigued  seismologists  for  many 
years  and  still  does.  As  a  first  and  subsequently  widely  accepted  substitute  for  EQ  energy  release, 
Richter  introduced  his  magnitude  scale  for  Californian  events  in  1935.  Over  the  years,  a  large 
number  of  magnitude  scedes  has  been  introduced  for  a  variety  of  seismic  phases  and  different 
tectonic  provinces  (e.g.,  see  Bath  1981,).  All  these  magnitude  scales  are  basically  similar,  namely 
of  the  form; 

^Vf  =  log(^)  +  !;(A,/)  +  c  (1) 

where  M  is  magnitude.  .4  is  peak  ground  displacement  for  the  seismic  phase  in  question,  y(A,/) 
incorporates  geometrical  spreading  and  inelastic  attenuation  terms  and  c  represents  the  station 
correction.  The  term  T  is  the  corresponding  signal  period. 

The  event  magnitude  parameter  is  coupled  to  the  seismic  moment  Mo  (Aki,  1967)  through  the 
formula  of  the  form; 


log  Wo  =  1.03W-  17.1 


(2) 


This  formula  is  taken  from  Sereno  et.  al.  (1988)  and  M  relates  to  Lg  phases  for  local/regional 
events.  Although  the  magnitude  parameter  is  less  clearly  related  to  the  earthquake  size  than  the 
moment,  the  magnitude  is  populzir  due  to  its  versatility  and  computational  ease.  Recently,  noise 
level  has  become  an  important  parameter  for  estimating  network  surveillance  capabilities,  say  in 
the  context  of  a  comprehensive  nuclear  test  ban  treaty.  In  this  study  we  focus  on  an  approach  to 
estimate  event  magnitude  in  near  real-time  and  start  with  a  brief  description  of  current  seismic 
signal  detectors.  Obviously,  parameters  extracted  during  this  operation  have  to  serve  as  a  basis 
for  magnitude  estimation. 

The  most  popular  signal  detectors  in  use  are  the  so-called  sliding  window  STA/LTA  type  which 
is  a  comparison  of  short  and  long  term  trace  amplitude  (a,)  averages  (e.g.,  Ruud  and  Husebye, 
1992).  Common  STA  definitions  are  of  the  forms; 


STA{abs) 


STA(rTns) 


±M 

tsl 


(3) 

(4) 


When  the  STA/LTA  ratio  exceeds  a  present  threshold,  the  presence  of  a  seismic  signal  is  declared. 
Our  preference  is  for  the  RMS  related  STA-LTA  definitions  since  signal  RMS  is  related  to  the 
signal  power  spectrum  via  the  Parseval’s  theorem,  namely; 


By  extensive  analysis  of  observational  data  Ringdal  and  Kvaerna  (1989)  have  established  empirical 
relationships  between  STA  and  event  magnitudes.  We  want  to  pursue  another  approach,  namely 
to  attempt  relating  the  STA(RMS)  to  maximum  ground  motion  amplitude  of  a  given  P/S/Lg 
-  phase  and  this  conform  with  the  original  Richter  (1935)  magnitude  definition.  This  problem 
is  rather  similar  to  a  one  dealt  with  in  seismic  htizard  analysis,  namely  how  to  generate  retdistic 


8 


% 


4ri 


•  • 


•  •  • 


•  • 


•  •  • 


ground  motion  spectra  using  random  vibration  theory  (RVT)  which  may  be  particular  appropriate 
if  the  waveforms  have  a  random  character  (Boore,  1983).  From  the  RVT  results  of  Cartwright 
and  Longuet-Higgins  (1956)  we  have 

M)  (6) 

where  ^(Aniaz)  is  the  expected  value  of  Amax  and  N  is  the  number  of  extremes  in  the  Arms  trace 
window.  When  N  is  large,  f{N)  can  be  approximated  by; 

/(Af)«(21niV)i  (7) 

When  N  becomes  small,  say  between  2  and  5,  a  better  approximation  may  be; 

/(N)«(21nlV)i +7(21n^V)-i  (8) 


where  7  is  Euler’s  constant  and  equals  to  0.5772 . . .  The  Eq{4)  and  Eq{6)  establish  a  relationship 
between  maximum  trace  amplitude  and  the  trace  RMS  value.  If  we  adopt  the  RMS  -variance  of  the 
STA-  definition  in  Eq(3),  then  an  estimate  of  maximum  signal/phase  amplitude  would  be  at  hand 
in  near  real-time.  For  the  conversion  to  event  magnitude  we  need  an  epicenter  distance  estimate 
which  in  case  of  arrays  may  be  tied  to  the  ’best’  beam  location.  Since  the  STA  parauneter  also 
is  available  for  the  non-detection  state,  that  is  for  the  pure  noise  portions  of  the  recordings,  the 
corresponding  maximum  amplitude  estimates  provide  a  measure  of  the  lower  magnitude  threshold 
for  seismic  event  detection.  We  can  also  determine  the  corresponding  error  in  the  magnitude 
estimate  by  a  simple  derivation  of  Eq{l); 


dM  = 


1  dA 
In(lO)  A 


(9) 


where  dA  is  the  derivative  of  the  A(moz)  value. 

In  RVT  theory,  derivation  of  Eqs(^  -  7),  the  time  series  is  presumed  stationary  and  individual 
samples  to  be  independent.  This  is  hardly  the  case  for  seismic  recordings,  so  we  have  simply  started 
to  test  validity  ot  Eqs(h -7)  using  real  seismic  recordings.  The  aim  of  such  sm  experiment  is  shown 
in  Fig.  2.2.1;  the  following  comments  apply.  There  is  clearly  a  simple  linear  relationship  between 
maximum  trace  amplitude  and  the  corresponding  RMS  value  for  the  entire  seismic  recordings. 
The  scattering  is  moderate  given  that  the  RMS  or  STA  window  length  is  at  least  twice  the 
dominant  signal  frequency.  In  case  of  band  filtered  records  (for  noise  suppression),  the  center 
frequency  of  the  pass  band  is  a  good  measure  of  the  dominant  signal  frequency.  The  number  of 
extremes  within  the  RMS  window  at  least  for  bandpass  filtered  records  is  also  very  stable.  The 
corresponding  magnitude  errors  according  to  Eq{%)  seldom  exceed  ±0.1  units  which  is  considered 
fully  3u:ceptable.  For  shorter  windows,  the  Amax ! ST A{RM S)  relation  still  holds  but  scattering 
increases.  We  conclude  in  the  basis  of  present  experience  with  analysis  of  real  recording  that  new 
real-time  event  magnitude  estimate  can  be  possible. 

We  are  now  undertaking  works  on  incorporating  the  ‘on-line  magnitude’  concept  in  the  auto¬ 
matic  analysis  of  seismic  network  recordings.  In  exrunple,  Ruud  and  Husebye  (1993)  have  recently 
demonstrated  that  the  3C  signal  detector  could  reliably  pick  phase  onset  times  while  correct  phase 
identification  could  occasionally  be  problematic.  We  expect  that  on-line  magnitude  estimates  for  a 
specific  event  location  would  contribute  towards  better  phase  identifications.  In  a  threshold  mon¬ 
itoring  context  we  could  have  an  attempt  to  on-line  estimation  of  seismic  network  event  detection 
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capability  using  on-line  noise  magnitude  estimates  from  globally  distributed  station  and  arrays 
(e.g.,  Sereno  et.  al.,  1988,  Kvserna  and  Ringdal,  1992). 
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Figure  2.2.1:  (a)  shows  a  sample  event  whose  frequency  bandwidth  is  between  0.75  -  1.5Hz. 
(b)  The  relation  between  the  RMS  (STA)  values  and  the  max  amplitudes  for  a  sliding  window 
is  illustrated  and  the  linear  line  which  fits  best  in  the  sense  of  residual  error  is  also  depicted. 
The  STA  window  length  is  3.2  sec  and  the  ssunpling  frequency  of  the  signal  equals  to  25Hz.  (c) 
shows  the  number  of  extremes,  and  (d)  the  corresponding  maximum  amplitude  ev2duations  are 
depicted;  i)  solid  line  is  the  original  maximum  amplitude,  ii)  dashed  line  represents  RVT,  whereas 
Hi)  dotted  line  is  the  evaluations  for  asymptotic  expressions.  The  magnitude  estimation  error  for 
(e)  the  RVT,  and  (f )  the  asymptotic  expressions  are  depicted. 
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2.3  Crustal  structure  of  Iberia 

The  regional  array  Sonseca  (ESLA)  near  Toledo,  Spain  is  now  operational,  and  in  this  context  we 
have  taken  the  opportunity  to  examine  the  crustal  structure  of  the  Iberian  Peninsula  in  terms  of 
Pn  and  Sn  velocity  distributions  as  derived  from  tomographic  inversions  of  local  bulletin  data.  A 
brief  description  of  this  work  is  given  below;  we  start  with  a  few  words  on  the  evolution  of  Iberia. 

The  geodynamic  evolution  of  the  Iberian  peninsula  and  the  western  Mediterranean  is  undoubt¬ 
edly  complex  (Dewey  et  al.,1989;  Platt  and  Vissers,  1989),  involving  both  compressional  tectonics 
related  to  the  convergence  of  the  African  and  Eurasian  plates  during  the  Eocene,  as  evidenced  by 
the  Pyrenees  and  the  Betic  Cordillera,  and  subsequent  extensional  tectonics,  as  evidenced  by  the 
creation  of  the  Valencia  Trough  and  the  Alboran  Sea.  It  is  expected  that  the  evolution  of  some  of 
these  major  tectonic  features  may  be  reflected  in  the  deeper  lithospheric  structure. 

Extensive  seismic  and  seismological  studies  of  the  Iberian  lithosphere  have  been  undertaken  over 
the  last  few  years.  Several  long-range  seismic  refraction  profiles  have  been  carried  out  traversing  the 
peninsula  and  at  sea,  providing  details  of  crustal  thickness  and  velocities,  both  at  its  Hercynian 
core  (e.g.  Medialdea  et  al.,  1986;  Cordoba  et  al,  1988;  Surinach  and  Vegas,  1988,  ILIHA  DSS 
Group, 1993)  and  in  areas  such  as  the  Betic  Cordillera  and  Alboran  Sea  which  have  special  tectonic 
significance  (e.g.  Medialdea  et  al.,1986,  Banda  et  al.,1993)(Fig.  2.3.1).  Additional  information 
can  also  be  obtained  using  data  collected  by  dense  seismograph  networks;  arrival  times  of  Pn  and 
Sn  phases  recorded  by  these  networks  can  be  antdysed  using  tomographic  techniques  to  obtsun 
estimates  of  sub-Moho  P-  and  S- velocities. 

In  this  study  we  derive  estimates  of  the  Pn  and  Sn  velocity  distribution  beneath  the  Iberian 
peninsula  by  carrying  out  tomographic  inversion  on  a  subset  of  earthquake  arrival  time  data 
routinely  collected  by  the  Spanish  National  Seismograph  Network  (Fig.  2.3.1b).  The  tomographic 
inversion  edgorithm  used  is  a  least  square  conjugate  gradient  (CGLS)  technique  described  in  detail 
by  Spakman  and  Nolet  (1988)  and  Blanco  and  Spakman  (1993).  In  the  inversion  we  allow  for 
known  variation  in  crustal  thickness  and  include  the  hypocentral  partuneters  as  unknowns.  This 
permits  further  event  relocation  using  the  velocity  estimates  found  in  the  inversion  (e.g.,  see 
Bannister  et  al.,  1991). 

We  use  a  two-layered  model  for  the  lithosphere,  the  upper  layer  of  which  is  taken  to  represent 
the  crust.  The  initial  model  of  the  crustal  thickness,  shown  in  Fig.  2.3.2.  was  primarily  based 
on  results  from  seismic  refraction  profiling  (e.g.  see  Badal  et  al.,1993)  supplemented  with  results 
derived  from  surface  wave  dispersion  studies  (e.g.  see  Badal  et  ai.,1992).  The  Iberian  crust  has  a 
relatively  uniform  thickness,  30-34  km,  but  is  notably  thicker  in  the  Pyrenees  (Surinach  et  al.,1993) 
and  the  Betics  (Banda  et  al.,1993).  The  crust  thins  rapidly  towards  the  Alboran  Sea;  Paulssen 
and  Visser  (1993)  give  a  thickness  of  only  20  km  in  the  coastal  area  around  Maletga. 

Forward  modelling  of  the  seismic  travel  times  was  initiedly  carried  out  using  ray-tracing  through 
a  starting  model  involving  no  latersil  velocity  variation,  only  crustad  thickness  variations.  Cricically 
refracted  Pn  and  Sn  phases  are  modelled  travelling  just  beneath  Moho  in  the  upper  mantle.  Both 
model  layers  are  sub-divided  into  blocks  of  constant  velocity  and  size  0.5°  x  0.5'  in  latitude  and 
longitude.  Ray  paths  along  great  circles  may  appear  slightly  curved  relative  to  the  chosen  grid 
system.  Variation  in  crustal  thickness  immediately  beneath  a  seismograph  station  or  event  was 
allowed  for  by  varying  the  thickness  of  the  upper  layer  (the  crust)  but  keeping  the  (local)  Moho 
interface  horizontal. 
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The  rays  were  traced  as  straight  lines  between  the  hypocentre  and  the  station  in  both  the  upper 
(crust)  and  lower  (sub-Moho)  layers.  Angies  of  incidence  at  the  layer  boundaries  were  taken  into 
account  using  the  respective  block  velocities.  The  effect  of  a  local  potentially  slanting  Moho  was 
ignored.  Sub-Moho  ray  paths  were  assumed  to  be  horizontal,  sJthough  differences  in  the  crustal 
thickness  beneath  the  station  and  hypocentre  were  accounted  for  by  introducing  a  length  scaling 
factor.  Initial  model  reference  P  and  S  velocities  were  6.4  and  3.65  km/s  in  the  crust  and  8.0 
and  4.57  km/s  sub-Moho  respectively.  These  Pn  and  Sn  velocities  were  derived  from  regression 
analysis  of  the  Pn  and  Sn  travel  time  observations  and  thus  reflect  the  velocity  model  used  in 
routine  hypocentre  determinations  for  the  Iberian  network. 

The  use  of  tomographic  row-action  techniques  prevents  us  from  cedculating  formal  resolution 
and  standard  errors  for  the  estimated  (unknown)  parameters.  As  a  substitute  we  inverted  a  syn¬ 
thetic  data  set  constructed  by  using  the  ray  population  from  the  actual  data  set  while  calculating 
travel  times  using  utiflcial  slowness  anomalies.  The  synthetic  model  used  consisted  of  alternating 
bodies  with  low  (8.0  km/s)  and  high  (8.40  km/s)  sub-Moho  velocity,  with  a  background  velocity 
of  8.2  km/s.  The  crust  was  presumed  to  be  homogeneous  and  with  an  initial  velocity  of  6.40  km/s. 
As  expected  the  velocity  anomalies  in  the  northern  parts  of  Spain  and  in  the  eastern  Pyrenees  are 
not  well  resolved  due  to  the  poor  ray-path  sampling.  In  contrast  the  P-velocity  reconstruction  for 
the  southern  and  western  parts  of  Iberia  appears  adequate,  as  it  does  also  for  parts  of  the  Ebro 
basin  and  the  Pyrenees.  The  reconstructed  velocity  pattern  smears  seaward  and  in  many  coastal 
areas,  except  to  the  west  of  Lisbonne  and  in  the  Gulf  of  Cadiz  where  it  appears  that  useful  results 
may  be  obtainable.  Similar  resolution  can  be  expected  in  the  S-wave  velocity  because  ray  paths 
are  similar. 

The  inversion  method  used  here  has  some  limitations,  as  noted  above.  One  assumption  is 
that  all  of  the  sub-Moho  arrivals  represent  head  waves,  with  no  regard  for  distance.  Synthetic 
modelling  (e.g.  Hestholm  et  al  '^4)  demonstrates  that  Pn  energy  may  represent  both  sub-Moho 
hesul  waves  and  turning  waves  within  the  upper  mantle.  It  is  difficult  to  differentiate  between 
these  two  wave  paths  using  only  travel  time  observations  and  there  is  little  knowledge  of  the  P- 
wave  velocity  gradients  within  the  upper  mantle  beneath  the  Iberian  peninsula.  The  head-wave 
assumption  should  be  reasonable  given  the  relatively  limited  distance  range  involved  for  most  of 
the  data.  It  is  difficult  to  get  a  quantitative  idea  of  the  degree  to  which  this  assumption  may  affect 
the  inversion  results. 

Fig.  2.3.3a  shows  the  Pn  velocity  distribution  determined  from  inversion  of  the  Pn  arrival 
time  data  set.  The  Pn  velocities  vary  between  7.7  km/s  (light  grey)  and  8.3  km/s  (dark  grey). 
The  striking  feature  is  the  low  (7. 7-7.8  km/s)  velocity  feature  in  the  Alboran  Sea,  to  the  west  of 
Gibraltar  Strait,  and  to  the  south  of  Portugal  in  the  Atlantic.  This  feature  is  presumably  related 
to  the  Neogene  extension  of  the  Alboran  Sea  region.  There  is  also  a  high  velocity  anomaly  beneath 
the  Sierra  Nevada  Ranges  (Betics).  The  sensitivity  modelling  shows  that  the  ray  path  coverage  in 
these  areas  is  adequate  for  good  resolution. 

The  Pn  velocity  beneath  the  Iberian  Massif  and  central  Spain  is  relatively  constant  at  around 
8.0-8. 1  km/s.  This  velocity  is  similar  to  that  found  in  refraction  studies  of  the  region  (e.g.  Banda 
et  al.,1981).  For  northwest  Spain  we  observe  slightly  lower  Pn  velocities,  around  7.9-8. 0  km/s, 
close  to  that  observed  from  several  seismic  refraction/wide-angle  profiles  carried  out  close  to  the 
Atlantic  coastline  in  that  region  (Cordoba  et  al.,1988). 

Shear  wave  (Sn)  velocities  estimated  from  the  tomographic  inversion  are  shown  in  Fig.  2.3.3b, 
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with  light  grey  shades  representing  slow  velocities  and  dark  grey  sh^ule8  rep'f  -^<‘uting  high  shear 
wave  velocities.  Low  Sn  velocities,  of  around  4.50  km/$,  are  observed  both  in  the  Albortui  Sea 
and,  onland,  beneath  Gibraltar,  as  well  as  beneath  Tagus  Basin  and  Central  System.  High  Sn 
velocities,  of  around  4.65  km/s,  are  observed  close  to  the  Guadalquiver  Basin  and  the  central  part 
of  the  Betic  Cordillera.  Little  velocity  variation  is  observed  across  most  of  the  Iberian  Massif. 
In  a  detailed  study  of  velocity  dispersion  of  surface  waves  ricraas  the  Iberian  peninsula  Badal  et 
al  (1992)  found  only  small  lateral  variations  in  shear  wave  velocities  over  most  of  the  peninsula. 
They  did  however  find  significant  effects  associated  with  the  toot  of  the  Betic  Cordillera,  with  low 
shear  wave  velocities  at  around  50  km  depth,  but  high  velocities  (relative  to  the  Iberian  Massif) 
near  the  Mobo.  It  is  clear  that  the  shear  wave  velocity  distribution  is  quite  complex  beneath  the 
Betic  Cordillera. 

In  summary,  Pn  and  Sn  velocity  anomalies  have  been  derived  for  the  Iberian  Peninsula  which 
correlate  with  some  of  the  tectonic  units  of  the  region.  Low  Pn  velocities  are  clearly  related  to 
the  Alboran  Sea  region  2uad  the  area  west  of  Gibraltar  Strait.  High  Pn  velocities  are  mainly 
related  to  parts  of  the  Hercinian  Massif  and  the  Betic  Cordillera.  Low  Sn  velocities  are  associated 
with  the  Alboran  Sea.  Our  results,  obtained  from  tomographic  inversion,  are  complementary  to 
results  obtained  from  a  variety  of  previous  studies  of  the  crust  and  mantle  beneath  the  Iberian 
peninsula  (e.g.  Mezcua  and  Udias,1991,  Badal  et  al.  1992,  Payoet  al.  1972  a,b,  Banda  et  Ed.,1993, 
Plomerova  et  al.  1993,  and  Blanco  and  Spakman,  1993). 


S.C.  BannUter  G.  Payo  E.S.  Husehye 

(Wellington)  (Toledo) 


References 

Aki,  K.  and  Lee,  W.H.K.,  1976.  Determination  of  three-dimensional  velocity  anomalies  under  a 
seismic  array  using  first  P  arrival  times  from  local  earthquakes  (part  1).  J.  Geophys.  Res..  81, 
4381-4399. 

Badal,  S.,  Corchete,  V.,  Payo,  G.,  Seron,  F.J.,  Canas,  J.A,  Pujades,  L.,  1992.  Deep  structure  of 
the  Iberian  Peninsula  determined  by  Rayleigh  wave  velocity  inversion,  Geophys.  J.  Int.,  108, 
71-88. 


Badal,  S.,  Gallart,  J.,  and  Paulssen,  H.,  (eds)  1993.  Seismic  studies  of  the  Iberian  Peninsula, 
Tectonophysics,  221,  1-34. 

Banda.E.,  Surinach.E.,  Aparicio.A.,  Sierra.J.,  and  Ruiz  de  la  Parte.E.,  1981.  Crust  and  upper 
mantle  structure  of  the  central  Iberian  Meseta  (Spain),  Geophys.J.R.astr.Soc.,  67,  779-789. 

Banda,  E.,  Gallart,  J.,  Garcia-Duenas,  V.,  Danobeitia,  J.J.  and  Makris,  J.,  1993.  Lateral  variation 
of  the  crust  in  the  Iberian  peninsula:  new  evidence  from  the  Betic  Cordillera,  Tectonophysics, 


14 


221,  53-60. 


Bannister, S.C.,  Ruud,B.O.,  and  Hiiaebye,E.S.,  1991.  Tomographic  estimates  of  sub-Moho  seismic 
velocities  in  Fennoscandia  and  structural  implications.  Tectonophysics,  189,  37-53. 

Bianco,  M.J.  suid  Spakman,  W.,  1993.  The  P-wave  velocity  structure  of  the  mantle  below  the 
Iberian  Peninsula:  evidence  for  subducted  lithosphere  below  southern  Spain,  Tectonophysics, 
221,13-34. 

Cordoba,  D.,  Banda,  E.  and  Ansorge,  J.,  1988.  P-wave  velocity-depth  distribution  in  the  Hercynian 
crust  of  NW  Spain,  Phys.  Earth  Planet.  Int.,  51,  235-248. 

Dewey,  J.F.,  Helmsin,  M.L.,  Turco,  E.,  Hotton,  D.W.H.  and  Knott,  D.S.,  1989.  Kinematics  of  the 
western  Mediterranean,  in  Alpine  Tectonics,  pp.  2654-283,  eds.  Coward,  M.P.,  Dietrich,  D. 
tmd  Park,  R.G.,  Geological  Society  of  London  Special  Publ.  45. 

Hestholm,  S.O.,  Husebye,  E.S.,  and  Ruud,  B.O,  1994.  Visualizing  seismic  wave  propagation  in 
complex  crust  -  upper  mantle  media  using  2D  finite  difference  synthetics.  Geophys.  J.  Int..  in 
pre». 

ILIHA  DSS  Group,  1993.  A  deep  seismic  sounding  investigation  of  lithospheric  heterogeneity  and 
anisotropy  beneath  the  Iberian  Peninsula,  Tectonophysics,  221,  35-51. 

Medialdea,  T.,  Surinach,  E.,  Vegas,  R.,  Banda,  E.  and  Ansorge,  J.,  1986.  Crustal  structure  under 
the  western  end  of  the  Betic  Cordillera  (Spain),  Ann.  Geophys.,  4(B4),  457-464. 

Mezcua,  J.  and  Udias,  A.,  1991.  Seismicity,  seismotectonics  and  seismic  risk  of  the  Ibero- 
Maghrebian  Region,  Monogratia  no.  8,  IGH,  Institute  Geogrsdico  Nacional,  Madrid,  Spain, 
3-15. 

Paulssen,  H.  and  Visser,  J.,  1993.  The  crustal  structure  in  Iberia  inferred  from  P-wave  coda, 
Tectonophysics,  221,  111-124. 

Payo,  G.,  1972a.  P-wave  residuab  at  some  Iberic  stations  and  deep  structure  of  south-western 
Europe,  Geophys.  J.R.  Astron.  Soc.,  26,  481-497. 

Payo,  G.,  1972b.  Crust-mantle  velocities  in  the  Iberitui  Peninsula  and  tectonic  implications  of  the 
seismicity  in  this  area,  Geophys.  J.R.  Astron.  Soc.,  30„  85-99 

Platt,  J.P.  and  Vissers,  R.L.M.,  1989.  Extensional  collapse  of  thickened  continental  lithosphere: 
a  working  hypothesis  for  the  Alboran  Sea  and  Gibraltar  arc,  Geology,  17,  540-543. 

Plomerova,  J.,  Payo,  G.  and  Babuska,  V.,  1993.  Teleseismic  P-residual  study  in  the  Iberian 
Peninsula,  Tectonophysics,  221,  1-12. 

Smith,W.H.F.,  and  P.Wessel,  1990.  Gridding  with  continuous  curvature  splines  in  tension.  Geo¬ 
physics,  55,  293-305. 


15 


•  • 


•  • 


•  • 


•  • 


•  • 


•  • 


•  • 


•  •  • 


Spakman,  W.  and  Noiet,  G.,  1988.  Imaging  algorithms,  accuracy  and  resolution  in  delay  time 
tomography.  In:  N.J.  Vlaar  et  al  (eds.),  Mathematical  Geophysics;  A  survey  of  recent  devel¬ 
opments  in  seismology  and  geodynamics,  Reidel,  Dordrecht,  155-188. 

Surinack,  E.  and  Vegas,  R.,  1988.  Lateral  inhomogeneities  of  the  Hercynian  crust  in  central  Spain, 
Phys.  Earth  Planet.  Int.,  51,  228-234. 

Surinack.  E.,  Marthelot,  J.-M.,  Gallart,  J.,  Daignieres,  M.  and  Him,  A.,  1993.  Seismic  images  and 
evolution  of  the  Iberian  crust  in  the  Pyrenees,  Tectonophysics,  221,  67-80. 

Weasel, P.,  rmd  W.H.F. Smith,  1991.  Free  software  helps  map  and  display  data,  EOS  Trans. AGU, 
72,  441,  445-446. 


•  • 


•  • 


•  •  • 


•  • 


•  • 


•  • 


•  • 


16 


Figure  2.3.1:  (a)  Tectonic  provinces  comprising  the  Iberi^ul  continental  area.  Symbols  are  the 
following:  1.  Deformed  Mesozoic  cover  of  the  Pyrenees,  2.  Mesozoic  area, ..  ilercynian  basement, 
4.  Mesozoic  external  units  of  the  Betics,  5.  Mesozoic  and  Paleozoic  internal  units  of  the  Betics,  6. 
Intermediate  Flysh  of  the  Gibraltar  arc,  7.  Undeformed  Mesozoic  cover.  8.  Tertiary  basins.  IM, 
Iberian  Massif.  P,  Pyrenees.  CCR,  Catalan  Coastsd  Range.  IC,  Iberian  Cordillera.  CS,  Central 
System.  B,  Betics.  DB,  Duero  Basin.  TB,  Tajo  Basin.  GB,  Guadalquivir  Basin.  EB,  Ebro  Basin. 
VT,  Valencia  Trough,  (b)  Seismograph  stations  that  reported  readings  used  in  the  analysis.  The 
Sonseca  (ESLA)  array  is  located  at  39.67“N,  3.96*W. 
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Figure  2.3.2:  A  model  of  Moho  depth  (in  km)  beneath  the  Iberian  peninsula,  gridded  using  GMT 
software  (Smith  and  Wessel,  1990). 
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Figure  2.3.3:  (a)  Pn  velocity  distribution  determined  from  inversion  of  the  observed  Pn  arrival 
time  data.  Dark  grey-scale  represents  high  velocities,  light  shades  represent  low  velocities,  in  km/s. 
(b)  Sn  velocity  distribution. 
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2.4  Crustal  modelling  for  2D  synthetics 


Structural  models  used  in  synthetic  seismogram  analysis  are  often  very  simple  in  order  to  avoid 
awkward  and  computationally  clumsy  mathematical  representations.  On  the  other  hsmd  complex 
models  could  be  generated  from  simple  ones  by  introducing  random  perturbations  in  physical  pa¬ 
rameters  like  velocity  and  density  or  rsindomizing  the  topography  of  the  free  surface.  Hestholm  et 
al  ( 1994)  among  others  have  used  von  Karman  functions  of  various  orders  for  simulating  hetero¬ 
geneous  lithosphere  models  for  2D  finite  difference  (FD)  synthetic  seismogram  experiments.  The 
results  here  are  encouraging  in  the  sense  that  models  incorporating  small  scale  velocity  perturba¬ 
tions  produce  realistic  seismograms  including  prominent  coda  wave  features  which  are  notoriously 
lacking  for  multilayered,  homogeneous  lithosphere  models.  In  general,  ray  synthetics  are  not  ob¬ 
viously  related  to  prominent  features  in  the  wavefield  recordings.  We  want  to  test  the  validity  of 
this  statement  by  computing  synthetics  for  such  multilayer  crustal  models  including  options  for 
introducing  perturbations  of  physical  parameters  as  mentioned  above.  In  other  context  we  want 
to  simulate  wave  propagation  across  prominent  tectonic  features  like  the  Viking  Graben  in  the 
North  Sea  which  acts  as  a  barrier  for  Lg-waves.  The  problem  in  common  here  is  to  generate  the 
proper  structural  models  in  a  simple  manner.  Our  approach  here  is  to  avoid  manual  digitalization 
to  the  extent  possible,  hence  our  preference  for  using  a  optical  scanner  in  combination  with  an  in¬ 
teractive  digitalization  system.  In  awidition  we  need  a  smart  algorithm  for  automatically  assigning 
grid  point  affiliations;  we  proceeded  in  the  following  manner. 

i.  Figure,  say  taken  from  a  publication,  is  prepared  for  the  scanner  like  the  MacIntoch  canvas 
system  with  output  (non-standard)  in  the  form  of  x,  y-coordinate  values. 

ii.  The  output  is  modified  in  such  a  way  that  the  model  'envelop'  haw  a  perfect  rectangular 
form  and  that  structural  boundaries  are  vectorized  being  represented  as  a  sequence  of  line 
segments.  Scaling  factors  are  incorporated  and  corrections  introduced  if  horizontal  and 
vertical  scalings  are  different. 

iii.  The  next  step  in  the  analysis  is  to  determine  all  closed  face  boundaury  curves  which  can  be 
established  from  all  boundary  line  segments  at  band  (step  ii).  For  example,  we  may  start 
in  the  lower  right  corner  in  Fig.  2.4.1  and  proceed  upwards:  at  the  first  line  intersection 
we  choose  the  leftmost  anti-clockwise  path,  and  so  on  until  ‘coming’  back  to  the  lower 
right  corner.  In  this  way  we  obtain  all  closed  face  boundary  curves  which  in  essence  are 
envelopes  for  areal  continuous  structural  elements  constituting  the  lithosphere  model.  Note, 
if  a  structural  element  is  completely  embedded  in  a  larger  one,  it  is  represented  by  two 
closed  curves  -  anti-clockwise  and  clockwise  respectively.  Also,  the  rectangle  encompassing 
the  whole  model  is  a  closed  curve  in  our  representation.  Note,  since  this  process  is  automatic 
there  is  no  formal  structural  element  identification,  except  for  a  single  1  to  N  numbering  of 
the  closed  curves  determined. 

iv.  The  last  and  most  tricky  step  was  to  introduce  both  a  numerical  element  identification 
and  in  puallel  tagging  such  a  label  to  all  the  respective  grid  points  contained  within  the 
proper  closed  boundary  curve.  Conventional  mathematical  geometry  proved  far  to  clumsy 
for  solving  this  problem,  but  recalling  a  theorem  from  integral  calculations  in  the  complex 
number  domain  a  simple  solution  was  found.  For  each  grid  point,  a  horizontal  branch  line  to 
the  left  was  formed.  Its  intersections  with  the  extremal  curves  (n,)  would  be  even  number 
for  a  point  outside  the  close  curve  and  an  even  number  for  a  point  inside.  The  identification 


number  (IDENT)  of  the  structural  element  to  which  the  grid  point  belongs  would  then  be 
a  sum  of  powers  with  base  number  of  2; 


N 

IDENT  =  (1) 

•=i 

where  N  is  the  number  of  closed  curves,  /(n.)  is  0  for  even  and  1  for  rti  odd.  Quite 
frankly,  we  were  amazed  over  the  very  simple  solution  to  the  above  rather  tricky  problem. 

V.  The  input  file  to  the  initial  digitalization  process  contained  a  listing  of  physical  parameters 
like  velocity  and  density  for  each  structural  element/layer  of  the  model.  All  this  information 
is  attached  to  a  single  point  within  each  individual  layer.  When  all  the  structural  elements 
have  been  identified,  it  is  easy  to  match  the  automatic  layer  labels  (eq.  1)  to  the  initial, 
deterministic  ones. 

The  workings  of  the  above  digitalization  procedure  is  demonstrated  in  Fig.  2.4.1.  An  added 
advantage  with  this  approach  is  that  input  file  errors  in  2D  PD  model  specifications  are  greatly 
reduced.  The  motivation  for  the  above  model  digitalization  venture  was  to  compute  synthetic 
seismograms  for  a  large  variety  of  lithosphere  models.  This  implies  that  the  subsequent  analysis  of 
synthetic  recordings  also  must  be  handled  in  an  efficient  manner.  Hence,  the  output  seismograms 
are  converted  to  SEG-Y  formats  and  then  subject  to  variou,:  kind  of  filtering,  stacking,  f-k,  p- 
tau-analysis  on  a  modern  workstation  like  ProMAX.  Also  this  task  has  been  accomplished  as 
demonstrated  in  Fig.  2.4.2. 

Future  plans  using  the  above  modelling  system  and  subsequent  waveform  analysis  tools,  in¬ 
volve  extended  testing  of  deep  crustal  models  derived  from  refraction  and/or  wideangle  reflection 
profiling  and  to  investigate  barrier  effects  like  the  apparent  blocking  of  the  Lg-propagation  across 
the  Viking  Graben  in  the  North  Sea. 
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Figure  2.4.1:  Model  of  P-wave  velocity  generated  by  the  above  procedure  and  used  for  2D  FD 
modelling.  An  explosion  source  was  placed  at  0.5  km  depth  at  x=10  km. 
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Figure  2.4.2;  ProMAX  screendump  of  synthetic  seismograms  generated  for  the  model  shown  in 
Fig.  2.4.1.  The  receivers  were  placed  with  3  km  intervals  in  the  range  x=13  to  x=280  km.  Time 
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3  Studies  completed 


In  this  section  we  gives  the  abstracts  for  three  papers  recently  completed  with  support  from  the 
present  AFOSR  Grant. 

3.1  Spatial  and  temporal  patterns  of  the  Fennoscandian  seismicity  -  an 
exercise  in  explosion  monitoring 

Matti  Tarvainen,  Institute  of  Seismology,  P.O.  Box  19,  FIN-00014  University  of  Helsinki.  Finland 

Eystein  S.  Husebye,  Institute  of  Solid  Earth  Physics,  University  of  Bergen,  Allegaten  41,  N-5007 
Bergen,  Norway 

Abstract 

Fennoscandian  stations  and  arrays  detect  daily  numerous  seismic  event  at  local  and  regionsd 
distances.  Almost  all  of  these  events  stem  from  various  kinds  of  explosions.  Such  explosion 
recordings  are  of  little  scientific  value,  but  add  considerably  to  the  daily  amalyst  workloads.  Also 
in  a  test  ban  verification  context  is  such  recordings  problematic  because  identification  between 
chemical  and  nuclear  explosions  is  not  always  obvious. 

In  this  work  we  report  on  efforts  to  automate  detection,  location  and  classification  of  seismic 
recordings,  stemming  from  local  mining  activities  and  other  explosions.  The  first  steps  here  are 
that  of  establishing  a  knowledge  base  on  mining  activities  such  as  their  locations,  firing  times, 
type  of  mining  activity  and  which  stations  are  most  likely  to  report  events  from  a  given  mining 
area.  Using  the  comprehensive  Helsinki  Bulletin  for  1991,  the  above  mentioned  parameters  have 
been  extracted  from  various  parts  of  Fennoscandia  and  north-western  Russia.  Consistent  diurnal, 
weekly  and  spatial  patterns  are  typical  for  most  of  the  seismic  events  included  in  the  bulletin. 

Comments;  Mining  and  other  kind  of  explosions  account  for  85  to  95%  of  the  entries  in  the 
seismic  bulletins  for  seismic  stations,  arrays  and  network  operations  in  Fennoscandia.  The  costs 
and  human  efforts  involved  in  proper  analysis  of  these  explosions  are  by  no  means  modest,  so 
the  question  arise  whether  this  observational  nuisance  possibly  could  be  handled  in  a  more  cost 
effective  manner.  For  example,  from  Fig.  2  in  the  above  study  (Appendix  A)  we  have  that  the 
explosion  activity  in  concentrated  to  a  few  areas  like  the  Estonian  oil  shale  quarries,  and  Karelian 
amd  Kola  mining  operations.  In  other  words,  we  try  through  cooperative  efforts  with  Finnish  and 
Russian  scientists  to  ensure  deployment  of  mobile  seismograph  equipments  for  a  period  of  a  few 
weeks  within  the  mentioned  mining  areas.  We  would  address  the  above  cost  efficiency  problem 
through  analysis  of  such  data  if  available  in  the  future. 


3.2  2D  finite  difference  elastic  wave  modeling  including  surface  topog¬ 
raphy 

Stig  0.  Hesthoim,  IBM  B«rgen  Environmental  Sciences  and  Solution  Centre,  Thormolensgate  55, 
N-5008  Bergen,  Norway 

Bent  O.  Ruud,  Institute  of  Solid  Earth  Physics,  University  of  Bergen,  Allegaten  41,  N-5007  Bergen, 
Norway 

Abstract 

A  2D  numericetl  finite  difference  silgorithm  accounting  for  surface  topography  is  presented.  Higher 
order,  dispersion-bounded,  cost-optimized  finite  difference  operators  ue  used  in  the  interior  of 
the  numerical  grid,  while  non-reflecting  absorbing  boundary  conditions  are  used  along  the  edges. 
A  transfer  from  a  curved  to  a  rect2Lngular  grid  achieves  the  modeling  of  the  surface  topography. 
Along  the  surface  we  use  the  free  surface  boundary  conditions.  To  account  for  the  surface  topog¬ 
raphy  is  important  for  the  complete  modeling  of  the  effects  of  wave  propagation.  Near-surface 
effects,  such  as  scattering,  are  otherwise  not  modeled  adequately.  Even  if  other  properties  of 
the  medium,  for  instance  randomization,  can  improve  numerical  simulations,  the  inclusion  of  the 
surface  topography  make  them  more  realistic. 

Comments:  Through  studies  of  teleseismic  events  we  have  earlier  found  that  scattering  by  surfatce 
topography  is  a  major  contributor  to  teleseismic  P-coda.  Also  for  regionad  and  local  events  surface 
waves  scattered  by  topography  can  be  found  in  the  low  frequency  part  of  the  P-  and  S-coda.  With 
the  present  method  we  a^e  able  to  model  such  effects  in  2D  with  little  extra  computational  toad 
compared  to  conventional  FD  methods.  An  extension  to  3D  is  possible,  but  presently  we  are 
prohibited  by  limited  computationad  capacity  to  do  realistic  3D  modelling. 
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3.3  Seismic  wave  propagation  in  complex  crust  -  upper  mantle  media 
using  2D  finite  difference  synthetics 
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Abstract 

Synthetic  seismogram  analysis  is  considered  an  adequate  tool  for  aiding  in  a  better  understand¬ 
ing  of  seismic  wave  propagation  in  the  lithosphere.  Our  finite  difference  (FD)  scheme  has  been 
used  for  computing  synthetic  seismograms  for  2D  crust/upper  mantle  models  of  size  150  x  400 
km~  and  with  options  for  a  free  surface  with  topography.  The  research  strategy  was  to  introduce 
successively  more  complex  lithosphere  models  for  generating  the  synthetics;  the  reference  model 
was  a  homogeneous  lithosphere.  The  extent  of  interface  scattering  is  visualized  through  displays 
of  synthetic  waveforms  and  snapshots  for  models  with  a  corrugated  Moho  only  and  free  surface 
topography  only.  The  latter  seems  to  dominate  and  in  the  form  of  P-to-Rg  and  S-to-Rg  conver¬ 
sions  at  the  surface.  Lithosphere  randomizations  were  introduced  through  von  Karman  functions 
of  order  0.3,  with  RMS  velocity  fluctuations  of  3-4  per  cent  and  correlation  distances  (horizontal 
and  vertical)  at  2.5  or  10  km.  In  case  of  a  medium  with  only  sub-.Moho  heterogeneities,  those 
with  horizontal  anisotropy  (ax  =  10  km;  az  =  2.5  km)  produced  relatively  strong  Pn  and  Sn 
phases.  The  respective  codas  were  dominated  as  in  most  of  our  experiments  by  P-to-S  and  S-to-S 
scattering  wavelets  excluding  Rg-scattering  at  a  free  surface  with  topography.  For  a  medium  with 
crustal  heterogeneities,  the  deformations  of  the  P  and  S  wavetrains  with  distance  were  clearly 
demonstrated.  For  full-scale  heterogeneous  lithosphere  models,  characteristics  features  of  the  syn¬ 
thetics  were  quantitatively  similar  to  observational  records  of  local  events.  Dominant  attributes 
were  a  pronounced  P  coda  consisting  mainly  of  P  and  Rg-scattered  wavelets,  and  a  relatively 
strong  S  coda  consisting  mainly  of  P-to-S  and  S-to-S  scattered  wavelets.  The  P  and  S  waveforms 
ttfe  severely  deformed  pointing  towards  the  futility  of  reliably  picking  many  secondary  arrivals  in 
local  event  recordings.  Most  of  the  scattering  wavelets  are  confined  to  the  crustal  waveguide  and 
to  surface  waves  since  coda  excitations  for  sensors  at  a  depth  of  100  km  were  weak  and  besides 
consisted  mainly  of  S  wavelets.  This  implies  that  a  strong  teleseismic  P  coda  does  not  reflect 
scattering  within  the  crust  in  the  source  region  but  preferably  a  complex  source  function.  Obser¬ 
vational  results  from  analysis  of  NORESS  and  ARCESS  bead  event  recordings  are  also  presented. 
Clearly  the  lithosphere  is  not  isotropically  inhomogeneous.  The  essence  of  our  2D  FD  synthetic 
seismogratm  experiments  is  that  a  simple  lithosphere  model,  being  moderately  heterogeneous,  gives 
rise  to  complex  seismograms  which  are  grossly  similar  to  the  observational  recordings.  In  contrast, 
complex  models  derived  from  profiling  surveys  (but  lawdeing  the  fine  scade  random  vairiations)  give 
simple,  ‘ray  tracing’  like  synthetics,  not  very  similar  to  observationad  records. 
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4  Future  work  plans 


The  focus  would  be  on  seismic  event  classification  and  network  monitoring  capabilities.  In  the 
former  case  we  would  also  explore  the  usefulness  of  broadband  3C  seismometer  recordings  and  to 
what  extent  SNR  enhancement  can  be  gauned  through  optimum,  adaptive  filtering  schemes.  A 
motivating  factor  here  is  that  inversion  of  low-frequency  waveform  data  provide  stable  estimates 
of  source  parameters  including  focal  depths. 

The  near  on-line  magnitude  estimation  scheme  may  prove  helpful  in  automating  seismic  network 
operations.  The  availability  of  amplitude  information  is  likely  to  improve  phase  identification  for 
non-array  stations  and  besides  almost  instantaneously  provide  estimates  of  test  ban  monitoring 
capabilities  on  local,  regionsd  and  global  levels. 

Most  of  our  knowledge  bearing  on  crustal  structures  and  compositions  stem  from  refraction 
and  wide-^ulgle  reflection  surveys.  A  basic  assumption  here  is  that  scattering  and  coda  excitation 
mechanisms  are  ignored.  We  are  planning  to  reexamine  published  results  derived  from  such  surveys 
in  Iberia  and  Kola  for  a  start  by  computing  2D  FD  synthetics  using  the  original  crustsd  models  in 
compuison  to  similar  models  being  subjected  to  random  perturbations  of  the  velocity  and  density 
parameters.  Likewise,  with  the  software  tools  at  hand  we  would  also  examine  the  structural 
peculiarities  acting  as  barriers  to  Lg-wave  propagation  say  across  the  North  Sea. 
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ABSTRACT 


Fennoscandian  stations  and  arrays  detect  daily  numerous  seismic  events  at  local 
and  regional  distances.  Almost  all  of  these  events  stem  from  various  kinds  of  explosions. 
Such  explosion  recordings  are  of  ^ttle  scientific  value,  but  add  considerably  to  the  daily 
analyst  workloads.  Also  in  a  test  ban  treaty  verification  context  is  such  recordings  prob¬ 
lematic  because  identification  between  chemical  and  nuclear  explosions  is  not  always 
obvious. 

In  this  work  we  report  on  efforts  to  automate  detection,  location  and  classification 
of  seismic  recordings,  stemming  from  local  mining  activities  and  other  explosions.  The 
first  steps  here  are  that  of  establishing  a  knowledge  base  on  mining  activities  such  as  their 
locations,  firing  times,  type  of  mining  activity  and  which  stations  are  most  likely  to  report 
events  from  a  given  mining  area.  Using  the  comprehensive  Helsinki  Bulletin  for  1991, 
the  above  mentioned  parameters  have  been  extracted  from  various  parts  of  Fennoscandia 
and  north-western  Russia.  Consistent  diurnal,  weekly  and  spatial  patterns  are  typical  for 
most  of  the  seismic  events  included  in  the  bulletin. 
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INTRODUCTION 


The  many  recordings  stemming  from  numerous  quarry  blasts,  mining  explosions 
and  other  man-made  seismic  events  are  a  nuisance  in  seismograph  network  operations  as 
their  analyses  add  considerably  to  the  daily  workload  of  the  analyst.  In  a  nuclear  test  ban 
context  an  additional  problem  is  that  event  classification  is  an  important  task  and  it  is  not 
always  simple  to  distinguish  between  earthquakes  and  explosions.  In  recent  years 
considerable  progress  has  taken  place  in  the  field.  Elaborate  schemes  based  on  so-called 
neural  network  techniques  have  been  proved  useful  (Dowla  et  al.  1990;  see  also  Tsvang 
et  al.  1993).  An  alternative  to  such  schemes  is  to  look  for  distinct  patterns  and  in  particu¬ 
lar  whether  signal  attributes  of  associated  seismic  recordings  are  site  specific.  For 
example,  Riviire-Barbier  and  Grant  ( 1991)  report  that  several  Fennoscandian  mines  have 
prominent  and  consistent  seismic  signatures,  although  variability  between  even  closely 
spaced  mines  (a  few  kilometers),  could  be  large.  The  idea  of  using  recording  signatures 
for  fast  reliable  event  location  and  source  type  identification  is  not  new  in  observational 
seismology  and  in  fact  is  daily  practiced  in  many  observatories  (e.g.  see  Vesanen  1944). 
Simply,  a  skilled  analyst  can  easily  differentiate  between  many  mining  areas  just  by  a 
quick  glance  on  the  relevant  wave  forms.  An  illustrative  example  here  is  that  the  Institute 
of  Seismology,  University  of  Helsinki,  more  than  two  decades  ago  introduced  a  template 
system  of  explosion  master  events  to  ease  the  daily  analysis  workload.  Essential  elements 
were  k  priori  knowledge  of  specific  mining  activities  and  which  explosions  produced 
distinct  and  easily  recognizable  recording  signals.  The  tabulation  of  this  coding  system, 
including  the  location  of  mining  operation  is  given  in  Table  1.  In  practice,  an  analyst 
wave  form  recognition  would  suffice  for  assuming  an  explosion  location  for  Che  event 
recordings  in  question.  There  is  now  a  desire  to  change  this  well-proven  but  essential 
manual  system  in  a  view  of  the  upgrading  of  the  national  Fennoscandian  networks  to 
digital  recording  and  the  general  need  for  operational  cost  efficiency. 

A  recent  approach  to  a  non-analyst  handling  of  the  mining  explosion  problem  is 
the  evolutionary  intelligent  monitoring  system  {IMS)  developed  for  automatic  multiar- 
rays  data  processing  analysis  (  Bache  et  al.  1990,  Bract  et  al.  1990).  Here  satellite  imagery 
information  has  been  incorporated,  providing  information  on  mine  locations  (including 
quarries)  and  other  types  of  large  scale  construction  works.  Such  information  has  proved 
very  useful  in  seismic  bulletin  work,  since  seismic  epicenter  locations  at  their  best  are  ac¬ 
curate  to  the  nearest  5-10  kilometers.  This  “spotting”  method  provides  accurate  epicenter 
solutions,  which  often  are  superior  to  these  derived  from  conventional  traveltime 
observations.  The  mining  and  quarry  sites  obtained  from  satellite  imagery  pictures  and 


those  of  the  old  manual  template  system  as  used  at  the  institute  of  Seismology,  Helsinki, 
are  shown  in  Figure  la  -  the  sites  coincide  well  with  each  other.  The  stations  in 
Fennoscandian  network,  which  contributed  data  to  this  study  is  shown  in  Figure  lb. 

In  this  paper  the  aim  is  to  track  the  Fennoscandian  mining  and  other  explosion  activity 
both  in  time  and  space  on  the  basis  of  available  bulletin  data  for  the  year  1991.  The 
distribution  of  those  events  is  shown  in  Fig.  2.  Also,  this  is  the  Hist  step  in  a  larger 
research  strategy,  to  discriminate  and  locate  automatically  and  accurately  regional 
seismic  events  from  digital  Fennoscandian  seismograph  stations. 

2.  OBSERVATIONAL  DATA  AND  DOMINANT  MINING  AREAS 

The  data  used  in  this  study  are  event  listings  from  Finnish  seismological  bulletins  and 
from  the  so-called  Nordic  Seismological  Data  Base.  The  latter  may  be  considered  a  re¬ 
gional  bulletin  as  these  event  listings  represent  a  merging  of  reports  from  the  seis- 
mological  observatories  in  Helsinki  (Finland),  NORSAR  (Kjeller,  Norway),  Bergen 
(Norway)  and  the  Russian  Academy  of  Sciences,  Kola  Branch  (Apatity).  The  Nordic  data 
base  at  present  is  seemingly  somewhat  incomplete  as  numerous  weak  events  detected  and 
automatically  located  by  the  NORESS  and  ARCESS  arrays  are  often  missing.  The 
epicenter  solutions  from  these  arrays  are  seemingly  not  officially  released  unless  checked 
by  the  local  analysts.  Rather  comprehensive  listings  for  ARCESS  and  NORESS  arrays 
are  given  by  Riviire-Barbier  ( 1993a,  1993b).  Anyway,  for  the  year  1991  altogether  4446 
event  reports  were  found  in  Helsinki  bulletins  (Uski  et  al.  1992),  within  the  geographical 
region  as  defined  in  Fig.  2.  Before  plotting  the  events  listings  were  prescreened  for 
removing  earthquakes  and  duplicate  events.  The  latter  are  defined  as  events  with  an 
origin  time  difference  less  than  30  s  and  epicenter  locations  within  50  km  of  each  other, 
but  different  reporting  agencies.  In  the  figure,  the  most  striking  feature  is  the  numerous 
explosions  in  the  Estonian  oil  shale  quarries,  from  which  1419  events  were  reported  in 
1991.  Other,  prominent  mining  areas  are  Kiruna-Malmberget  (northern  Sweden),  Russian 
Karelia,  and  the  Kola  district  (NW  Russia).  Explosion  activity  at  west-coast  of  Norway 
(Bergen)  is  tied  to  road  and  construction  works,  and  hence  has  a  more  dispersed  epicenter 
panem.  The  prominent  explosion  areas  were  subject  to  more  detailed  analysis,  presented 
below. 


2.1.  KmUNA-MALMBERGET  NORTHERN  SIVEDEN  (Figures  3  -  6) 


The  high  columns  in  Figure  3  coincide  with  the  Maimberget  mine.  There  is  also 
high  mining  activity  in  Kiruna.  The  explosion  activity  is  relatively  high  and  the  shooting 
times  are  mainly  in  the  afternoon  and  late  evenings  (Fig.  4).  These  are  two  distinct 
shooting  intervals  (Fig.  4),  which  may  reflect  different  shooting  practice  at  different 
mines.  To  see,  if  we  could  resolve  this  problem,  the  explosion  locations  were  plotted  for 
specific  time  intervals  (Fig.  S).  The  majority  of  events  occurring  in  afternoons  can  be 
connected  to  the  Maimberget  and  Aitik  mines  while  very  early  morning  events  are 
evidently  explosions  in  the  Kiruna  mine. 

The  Kiruna-Malmberget  mining  operations  are  unique  in  the  sense  that  reportings 
from  the  Swedish  seismograph  stations  are  lacking  as  shown  in  Fig.  6.  This  does  not 
reflect  poor  station  performances  but  instead  a  general  lack  of  economical  resources  for 
daily  analyst  analysis  of  local  (mining)  seismic  events.  The  mentioned  non-reportings  of 
the  NORESS  and  ARCESS  arrays  are  also  obvious  from  the  figure.  In  other  words,  the 
local  event  reporting  practice  between  Nordic  seismological  observatories  still  vary 
considerably.  As  a  rule  of  thumb  the  stations  closest  to  seismic  source  areas  should 
report  many  more  events  than  far  away  stations,  say  along  the  west  coast  of  Norway 
where  besides  the  background  noise  is  relatively  high.  So  far,  we  have  not  tried  to 
establish  sort  of  seismic  recording  attributes  (Joswig  and  Schulte-Theis  1993) 
differentiating  between  the  various  mining  operations  in  the  Kiruna-Malmberget,  because 
the  incomplete  bulletin  data  are  not  convenient  for  doing  so.  However,  the  importance  of 
this  issue  is  obvious  from  Fig.  5.;  the  conventional  epicenter  locations  are  not  always 
close  to  the  known  mining  locations. 

2.2  THE  RUSSIAN  KARELIA  AND  ADJACENT  AREAS  (Figures  7  -10) 

The  mining  activity  here  is  large  as  demonstrated  in  Figure  7  and  8.  where  also 
Estonian  oil  shale  quarry  explosions  are  included.  The  seismic  activity  has  been  subject 
to  special  studies,  e.g.  Tarvainen  (1992)  concluded  that  on  the  basis  of  event  locations 
tied  to  nearby  3-component  station  recordings  it  was  impossible  to  separate  individual 
mining  operations  (see  also  Bache  et  al.  1990).  Mines  identified  via  SPOT  satellite 
images  were  found  to  be  close  to  the  epicenters  (Fig.  8). 

The  diurnal  shooting  pattern  for  this  area  is  rather  distinct  with  a  strong  concentration  in 
the  afternoon  hours  on  working  days  (Fig.  9).  Explosion  charges  appear  to  be  small  since 
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the  rapottiag  frequency  map  in  Figure  10  is  dominated  by  Finnish  stations  tHeisinki 
agency).  A  few  reports  stem  from  Norwegian  arrays  and  Russian  stations. 
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23.  KOLA  PENINSULA  AND  ADJACENT  AREAS  (FigBrca  11-13) 


The  seismic  event  locations  reflect  the  diverse  and  extensive  mining  activities 
taking  place  in  this  area  (Fig.  1 1).  Again,  there  is  a  good  diurnal  clustering  of  events  in 
the  afternoon  and  there  is  an  increased  amount  of  explosions  in  weekends  (Fig.  12). 

Recordings  of  the  Kola  events  stem  mainly  from  the  Finnish  stations  and  the 
Helsinki  Data  Center  (Fig.  13).  More  extensive  studies  of  the  Kola  Peninsula  seismicity 
are  presented  by  Kremenetskaya  (1991),  Kremenetskaya  and  Tijapitsin  (1992),  and 
lately  by  Mykkeltveit  ( 1993). 

2.4.  WESTERN  NORWAY  EXPLOSIONS  (Figures  14-16) 

There  is  no  mining  operation  in  the  area  shown  in  Figure  14,  but  road  and 
construction  related  explosions  are  numerous.  192  events  occurred  around  Bergen  out  of 
which  46  were  listed  as  earthquakes.  This  explains  partly  the  dispersion  of  epicenters,  but 
also,  because  there  is  no  analyzed  template  system  in  use  at  the  Bergen  Observatory.  The 
diurnal  and  weekly  distribution  of  origin  times  is  also  relatively  dispersed  (Fig.  IS), 
eventhough  clear  concentration  at  IS  hours  on  working  days  is  clear.  Most  of  the  noc¬ 
turnal  events  are  earthquakes.  The  main  part  of  detecting  stations  is  located  in  western 
Norway  (Fig.  16). 

3.  DISCUSSION 

The  seismicity  in  Fennoscandia  is  very  moderate  with  two  or  three  earthquakes  a 
week,  equivalent  to  a  few  percents  of  our  explosion  data  base  for  1991.  Even  if 
teleseismic  events  are  taken  into  account,  the  dominance  of  local  explosions  prevails, 
contributing  80-9S%  of  the  total  event  recording  population  at  every  station  in  this 
region.  If  all  event  records  were  subject  to  the  same  careful  screening,  about  90  per  cent 
of  the  analyst  workload  would  be  tied  to  local  explosion  recordings  of  little  or  no 
scientific  value.  In  practice  such  work  allocations  are  not  feasible  over  extended  periods 
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of  times  at  most  seismological  observatories.  However,  as  most  of  the  explosion 
epicenters  are  located  in  nonhem  and  eastern  Fennoscandia,  the  excellent  bulletin  work 
performed  at  the  Helsinki  Data  Center  has  in  ^  '^ctice  proved  to  be  most  beneficial  for 
similar  work  at  the  other  Nordic  seismological  centers.  All  event  records  from  the  stations 
in  the  coastal  areas  of  Norway  (Bergen  Observatory)  are  carefully  analyzed,  but  the 
number  of  recorded  events  is  low  due  to  a  combination  of  high  noise  levels  and  also, 
relatively  high  signal  attenuation  for  travel  paths  beneath  the  (Caledonian  mountain  range. 
NORESS  and  ARCESS  produced  automatic  bulletins  in  1991,  but  their  event  listings 
were,  as  mentioned,  only  occasionally  included  in  the  Nordic  bulletins. 

The  major  outcome  of  our  study  of  the  Fennoscandian  seismicity  for  1991  is  that 
the  event  population  is  completely  dominated  by  numerous  mining,  quarry  and  other 
explosions.  Such  industrial  operations  are  often  stationary  in  space,  and  as  demonstrated 
exhibit  well  defined  temporal  patterns.  The  former  feature  is  not  always  obvious  from  the 
bulletin  data  as  illustrated  in  Fig.  5  for  the  Kiruna-Malmberget  mining  district  Note,  that 
most  of  these  event  locations,  do  not  stem  from  the  national  Finnish  bulletins,  but  from 
readings  of  stations  in  northern  Norway.  Anyway,  the  issues  of  incorporating  well- 
established  seismicity  information  and  stationary  signal  attributes  are  of  particular  im¬ 
portance  for  automatic  seismic  network  operations.  The  observational  data  at  hand  are 
often  less  voluminous  and  less  reliable  than  those  extracted  through  analyst  inspections  of 
station  records  (Ruud  and  Husebye  1992,  Ruud  et  al.  1993).  As  mentioned  above,  formal 
epicenter  solutions  could  be  unreliable  when  based  on  P-arrival  times  only  due  to  odd 
network  configurations.  The  kind  of  seismicity  information  to  be  incorporated  will 
naturally  be  spatial  and  temporal  characteristics  of  mining  operations.  For  example,  being 
confident  that  given  event  recordings  stem  from  a  specific  mine,  the  coordinates  of  that 
mine  should  replace  the  formal,  seismic  epicenter  solution. 

Another  seismicity  parameter  of  interest  is  station  detectability,  that  is  which 
stations  are  most  likely  to  record  explosions  from  a  specific  mine.  Unfortunately,  the 
bulletin  listings  are  as  mentioned  incomplete  on  this  account,  so  we  only  plotted  the  event 
reporting  frequency  for  those  stations  and  agencies,  for  which  data  were  available  (Figs. 
6,  10,  13  and  16).  Notice  that  the  dominance  of  reporungs  from  Helsinki,  reflects  the  ex¬ 
tensive  usage  of  the  template  system  in  Finland  and  does  not  refer  to  a  particular  station. 
This  kind  of  information  would  be  useful  for  detecting  phase  association  errors;  for 
example  western  Norway  stations  are  unlikely  to  detect  mining  explosions  in  the  Kola 
and  Karelian  areas.  Moreover,  most  of  the  Fennoscandian  mining  explosions  are  small 
magnitude  (Ml<2)  events  and  seldom  recorded  beyond  distances  of  a  d-5(X)  km  (Figs. 
6,10,13  and  16;  Ruud  and  Husebye  1992). 


To  use  signal  attributes  in  automatic  seismic  network  operations  require  waveform 
stationarity  for  a  given  station  for  a  specific  mine.  This  appears  often  to  be  the  case  in 
Fennoscandia  so  the  challenge  here  is  to  mimic  analyst  template  operation  numerically. 
The  mentioned  elaborated  IMS-system  (Bache  et  al.  1990,  Bratt  et  al.  1990,  Bache  et  al. 
1993)  has  provisions  for  implementing  such  features  and  also  learning  abilities  useful 
solutions  to  the  above  kind  of  “artiftcial  intelligence”  (Al)  problems  have  proved  difficult 
in  practice.  Anyway,  basic  requirements  here  are  compilation  of  reference  events  for 
every  station  in  the  network,  and  that  the  initial  of  epicenter  location  is  accurate 
(Rividre-Barbier  and  Grant  1993).  Besides,  no  network  can  be  operated  fully  automatic  in 
terms  of  producing  high-quality  bulletins,  simply  because  signal  attributes  for  weak 
events  would  be  unreliable  due  to  noise  contaminadon. 


4.  SUMMARY  REMARKS 


We  have  in  this  study  examined  the  seismic  aspect  of  the  mining  activities  and 
other  in  Fennoscandia,  which  are  characterized  by  clear  temporal  and  spatial  patterns. 
Analysts  in  Helsinki  have  taken  full  advantage  of  such  knowledge  so  as  to  significantly 
ease  the  daily  work  load. 

With  digital  seismograph  networks  being  deployed  in  various  parts  of  the  region 
the  challenge  is  to  teach  our  computers  to  locate  reliably,  and  identify  the  many 
thousands  of  explosions  being  recorded  annually.  This  would  require  much  painstaking 
work,  since  surprisingly  little  detailed  information  is  available  on  many  mining 
operations  in  particular  in  north-western  Russia,  but  also  elsewhere  (short  term 
construction  works).  However,  we  feel  confident  that  this  kind  of  problem  can  be  solved 
using  robust  “grid  search”  techniques  for  incorporating  i  priori  seismicity  knowledge  of 
the  kinds  dealt  with  here.  The  difficult  problem  is  to  mimic  numerically  the  analyst’s 
ability  to  quickly  recognize  mine  specific  waveforms.  This  implies  that  signal  attributes 
are  spatially  stationary,  which  is  often  the  case(Rivi^re-Barbier  and  Grant,  1993). 

Observationally  seismology  has  progressed  tremendously  over  the  last  decade  in 
terms  of  station  design  and  deployment,  digital  recording,  satellite  communication  and 
advanced  data  center  facilities  in  various  countries.  A  paradox  in  this  context  is  the  scant 
attention  paid  to  the  pressing  problem  of  automatic  signal  analysis  and  bulletin  pro¬ 
duction.  With  a  few  exceptions  too  small  efforts  are  invested  in  this  most  essential  aspect 
of  network  operation  and  earthquake  monitoring. 
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Table  1 


The  Helsinki  analyst  template  system;  latitude  and  longitude  of  the  operation.  If  known, 
the  corresponding  name  of  the  mine  is  given. 


Mines  and  explosion  sites  in  Estonia  and  Russian  Karelia 


Latitude 

mjjuj 

Name  of  mine 

Latitude 

Name  of  mine 

59.24 

— 

24J3 

59.41 

24.59 

59_5 

25.0 

5933 

27.27 

5933 

27.07 

59.27 

27.73 

59.24 

27.83 

5931 

27.63 

59.3 

27.5 

5937 

28.53 

59.45 

26.49 

60.0 

29.9 

60.02 

29.74 

59.6 

30.0 

5937 

28.43 

5936 

28.37 

61.01 

29.04 

60.90 

2935 

60.8 

29.3 

60.95 

29.18 

61.1 

30.3 

61.5 

30.4 

61.4 

31.6 

60.7 

29.0 

61.4 

343 

61.9 

30.6 

61.14 

29.87 

62.2 

343 
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64.68 

30.66 

RHIHHi; 

60.6 

29.2 

60.7 

28.7 

60.8 

29.5 

Norwegian  mines 


Latitude 

Name  of  mine 

Latitude 

Lonjtitude 

Name  of  mine 

5830 

6.40 

Titania 

69.6 

29.9 

Bjemavattnet 

Mining  sites  on  Kola  Peninsula  and  adjacent  areas 

Latitude 

Name  of  mine 

Latitude 

Name  of  mine 

67.67 

33.74 

Kirovsk 

67.63 

33.84 

Rasvumschorr 

67.64 

34.02 

Koashva 

67.64 

33.88 

69.4 

33.18 

68.16 

33.18 

67.6 

34.2 

69.6 

32.2 

67.7 

31.4 

67.56 

30.44 

Kovdor 

69.3 

34.4 

69.2 

34.7 

68.87 

33.03 

Murmansk 

69.23 

33.17 

Mines  in  Northern  Sweden 

Latitude 

Name  of  mine 

Latitude 

Name  of  mine 

67.18 


20.67 


67.12 


20.90 


Aitik 


67.8 


20.2 


Kiruna 


67.7 


21.0 


Svappavaara 

Mines  and  explosions  sites  in  Finland 


Latitude 

Loaxitude 

Name  of  mine 

Ladtude 

Name  of  mine 

60.17 

23.84 

Mustio 

6133 

23.03 

Stormi 

6030 

22.29 

Parainen 

61.9 

21.5 

Otamo 

61.6 

21.7 

Kridsiceh 

62.07 

27.41 

Ankele 

63.12 

27.72 

64.12 

28.06 

64.1 

24.7 

Perkkdd 

62.83 

29.25 

Horsmanaho 

63.66 

26.05 

62.8 

22.9 

Tomava 

63.16 

27.99 

Nilsia 

63.0 

26.8 

Talluskanava 

55.78 

24.70 

63,85 

25.05 

Hitura 

61.64 

24.26 

Orivesi 

61.94 

29.03 

Ruokojarv'i 

Fig.  la  .The  mining  sit 
at  the  Institute  of  Seisi 
system  (IMS)  based  < 
locations  overlap  nice! 


Fig.  2.  Seismic  events  of  presumed  explosion  origin  in  Fennoscandia  in  1991  according 
to  the  Nordic  Bulletin  issued  by  the  Finnish  National  Data  Center  (FNDC).  The  size  of  a 
column  is  1”  x  1”.  The  highest  activity  is  on  the  Estonian  coast,  where  many  oil  shale 
quarries  are  located.  The  highest  column  in  the  area  represents  145  events.  The  total 
number  of  explosions  for  the  Estonian  mining  operations  was  1419  in  1991.  Also,  in 
northern  Sweden,  western  Norway,  Russian  Karelia  and  Kola  the  explosion  activities  are 
high.  Most  events  in  Finland  stem  from  mines,  quarries  and  construction  work  sites. 
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Fig.  3.  Events  in  Kiruna-Malmberget  region  northern  Sweden,  in  1991.  The  mine  sites  are 
shown  by  black  dots.  The  black  diamonds  signify  small  earthquakes.  The  column  size  is 
0.  l*x  0. 1®  and  the  highest  column  represents  26  events.  The  total  number  of  events  in  the 
area  was  212  including  1 1  earthquakes. 


Fig.  4.  Daily  and  weekly  distribution  of  seismic  events  in  northeni  Sweden.  The  majority 
of  events  occun  on  working  days  early  in  the  morning  (local  time).  Daytime  saving  time 
from  spring  to  early  autumn  is  taken  into  account.  Afternoon  is  the  time  of  mining 
explosions  in  Aitik  and  Malmberget.  The  Kiruna  mining  explosions  occur  daily  (peak  on 
Saturday  mornings).  The  highest  column  is  equivalent  to  22  events. 


18  19  20  21 

Longitude  (East) 


22  18  19  20  21 

Longitude  (East) 


Fig.  S.  Distribution  of  e;tplosioas  in  northern  Sweden  for  different  diurnal  time  intervals . 
The  day  is  divided  into  four  intervals  of  six  hours.  The  shooting  practice  in  the  Kiruna 
mine  (Fig  5a)  is  clearly  different  from  that  in  the  other  mines  of  the  area. 


Fig.  6.  Event  reporting  frequencies  for  agencies  and  stations  in  Fennoscandia  and  Russia 
of  explosions  in  northern  Sweden.  The  Helsinki  Data  Center  and  stations  in  northern 
Norway  had  the  highest  repotting  rates  of  the  events,  but  without  identifying  individual 
detecting  stations.  The  mine  locations  are  plotted  as  black  dots. 


Fig.  7.  Events  in  Russian  Karelia  in  1991  and  adjacent  areas.  The  column  size  is 
0.1'*x0.1'*.  The  sites  of  highest  activities  match  well  with  mining  locations  seen  on 
satellite  imagery  data  (GSE/US/73).  The  outstanding  feature  is  the  many  explosions  in 
the  oil  shale  mines  in  eastern  Estonia,  which  for  purpose  of  comparison  were  included. 
The  highest  column  here  represents  51  events.  The  number  of  events  within  the  mapping 
area  was  488. 
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Fig.  8.  The  events  of  Fig.  7  shown  together  with  SPOT  located  mining  sites. 
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Fig.  9.  Diurnal  and  w< 
adjacent  areas.  They  ar 
during  prime  working  b< 


Fig.  10.  Event  reporting  frequencies  for  agencies/stations  of  explosions  in  Russian 
Karelia.  The  Finnish  stations,  most  of  which  are  within  300  km  exhibited  the  best 
monitoring  capabilities.  The  column  height  at  the  Helsinki  Data  Center  represents  416 
events,  which  is  close  to  86%  of  all  the  events  in  the  area  in  1991. 


Fig.  1 1.  Events  on  Kola  Peninsula  and  adjacent  areas.  See  fig.  3  for  symbols.  The  highest 
concentrations  occurred  near  known  mining  sites  (black  dots).  The  total  number  of  events 
in  the  area  was  780  for  the  year  1991, 1  he  highest  column  represents  38  events.  The  few 
occurring  earthquakes  are  marked  by  black  diamonds. 


Fig.  14.  Western  Norway  events  in  1991,  see  Fig.  3  for  symbols.  The  most  active  zone 
was  around  Bergen,  (the  highest  column  represents  only  7  events).  There  different  kinds 
of  construction  works  took  place.  The  area  contributed  192  events  in  1991,  out  of  which 
46  probably  were  listed  as  earthquakes. 


Fig.  15.  Diurnal  and  weekly  distribution  of  seismic  events  in  western  Norway.  The 
highest  activity  is  during  working  hours.  Typically,  there  is  less  activity  during  weekends. 
The  relatively  large  number  of  ea.  >  quakes  (46  out  of  191  events)  explains  the  dispersion 
in  the  diurnal  distribution. 
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Fig.  16.  Event  reporting  frequencies  for  agencies  and  stations  of  explosions  and 
earthquakes  in  western  Norway.  The  western  Norwegian  network  reported  mainly 
Norwegian  coastal  events,  with  the  highest  column  for  the  station  ASK  accounting  for 
241  events.  Some  events,  reported  also  in  Finland  and  western  Russia,  were  evidently 
earthquakes. 


A-Zl 


•  • 
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Abstract 

A  2D  numerical  finite  difference  algorithm  accounting  for  surface  topography  is 
presented.  Higher  order,  dispersion-bounded,  cost-optimized  finite  difference  op¬ 
erators  are  used  in  the  interior  of  the  numerical  grid,  while  non-reflecting  absorbing 
boundary  conditions  are  used  along  the  edges.  A  transfer  from  a  curved  to  a  rect- 
amgular  grid  achieves  the  modeling  of  the  surfture  topography.  Along  the  surface 
we  use  the  free  surface  boundary  conditions.  To  account  for  the  surface  topog¬ 
raphy  is  important  for  the  complete  modeling  of  the  effects  of  wave  propagation. 
Neair-surface  effects,  such  as  scattering,  are  otherwise  not  modeled  adequately. 
Even  if  other  properties  of  the  medium,  for  inst^mce  randomization,  can  improve 
numerical  simulations,  the  inclusion  of  the  surface  topography  make  them  more 
readistic. 


Introduction 

Including  topography  at  the  free  surface  of  an  elastic  medium  leads  to  improved 
modeling  of  near-surface  effects,  especially  those  in  the  high  frequency  part  of 
the  wave  field.  Experimentation  with  other  medium  characteristics  can  partly 


contribute  to  the  same  effects.  Nevertheless,  without  including  the  topography  in 
the  medium  model,  the  results  of  these  experiments  may  be  unreliable  with  regard 
to  near-surface  effects. 


In  recent  years,  finite  difference  methods  have  gained  popularity  among  geo¬ 
physicists  as  a  flexible  tool  for  modeling  seismic  wave  propagation  in  various  ge¬ 
ological  environments.  Although  most  heterogeneous  media  can  be  effectively 
handled  this  way,  the  inclusion  of  topography  at  the  surface  of  an  elastic  medium 
has  proven  to  be  a  non-trivial  problem  in  the  synthetic  calculations.  Finite  dif¬ 
ference  and  finite  element  methods  of  lower  order  have  often  been  used  for  the 
modeling  of  waves  in  the  vicinity  of  surface  topography.  The  diflBculty  with  the 
finite  difference  modeling  lies  in  adapting  the  grid  required  by  the  method  on  to 
the  topography  at  the  free  surface.  The  finite  element  method  is  the  more  flexible 
in  incorporating  variable  geometry,  but  it  is  more  time  consuming  for  the  same 
order  (Lysmer  L  Drake,  1972). 

A  cost-effective  method  well  suited  for  handling  different  types  of  geometries, 
is  the  boundary  element  method  e.  g.  see  Mansur  ic  Brebbia,  1981,  1982.  This 
method  applied  in  the  elastic  case  can  be  found  in  Cruse  k  Rizzo  (1968),  Cruse 
(1968)  and  Hall  k  Robertson  (1989)  among  others.  For  relatively  simple  structures 
in  the  interior  of  the  medium  (homogeneous  or  linearly  varying),  this  method  is 
advantageous  in  demanding  less  computer  storage  than  other  methods.  The  reason 
is  that  calculations  are  only  kept  track  of  along  the  boundaries.  The  boundaries  cem 
also  have  complex  shapes,  as  long  as  they  consist  of  closed  curves.  Nevertheless, 
full  matrices  have  to  be  inverted. 

A  procedure  for  modeling  free  surface  topography  with  lower  order  finite  dif¬ 
ferences  was  given  by  Jih,  McLaughlin  k  Der  (1988).  This  work  classified  different 
types  of  slopes  and  abrupt  changes  of  slopes,  and  prescribed  a  procedure  for  h2ui- 
dling  each  of  them  by  implementing  the  free  surface  conditions.  To  the  contrary 
of  Jih  et.  al.  (1988),  Tessmer,  Kosloff  k  Behle  (1992)  assume  a  smooth  surface, 
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i.  e.  the  first  derivative  of  the  topography  function  exists.  A  stretched  (or  curved) 
coordinate  system  is  initi2dly  applied,  and  its  upper  boundary  coincides  with  the  to¬ 
pography.  A  transformation  onto  a  rectangular  coordinate  system  is  done  in  order 
to  perform  the  computations,  using  the  Chebyshev  method  and  the  Fourier  method 
for  spatial  differentiation.  The  present  work  applies  this  transformation  from  a 
curved  to  a  rectangular  system,  while  using  higher  order,  dispersion-bounded, 
cost-optimized  finite  difference  operators  (Sguazzero,  Kindelan  k.  Kamel,  1989) 
for  discretizing  the  spatial  derivatives.  An  additional  coordinate  transformation 
has  to  be  done  locally  at  each  point  on  the  surfaM:e,  so  that  the  vertical  coordinate 
coincides  with  the  normal  vector  before  the  implementation  of  the  free  surface 
boundary  conditions.  The  discretization  of  these  are  done  by  lower  order,  centrad, 
staggered  finite  differences  (Fornberg,  1988,  a). 

Although  the  method  of  Tessmer  et.  al.  (1992)  for  incorporating  topography  at 
the  free  surface  cannot  be  used  directly  in  our  context  of  finite  difference  methods, 
it  proved  instructive  for  our  modification  of  the  boundary  conditions.  The  present 
work  concentrates  on  the  theory  behind  our  procedure  for  modeling  elastic  waves. 
In  forthcoming  works  results  and  analyses  will  be  presented,  as  well  as  comparisons 
with  a  long  range  of  real  measurements. 

The  following  sections  present  the  basic  equations  together  with  their  dis¬ 
cretization  onto  a  staggered  grid  (Sguazzero  et.  al.,  1989,  Levander,  1988).  The 
equations  for  the  curved  system  transferred  onto  the  recttmgultur  grid,  where  the 
computations  are  done,  will  be  presented.  Then  the  free  boundary  conditions  re¬ 
sulting  from  the  rotated  coordinate  systems  at  the  surface  are  given.  Next,  we 
indicate  how  the  numerical  discretization  is  achieved.  Finally,  we  present  a  few 
numerical  examples  with  a  random  topography  and  realistic  medium  P-velocity 
models.  These  are  illustrative  of  the  importance  of  near-surface  effects  for  the 
high-frequency  part  of  the  wave  field. 


•  •  • 
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Elastic  Wave  Modeling  Formulation 

The  basic  equations  governing  wave  propagation  in  a  continuous  elastic  medium 
ue  the  momentum  conservation  and  the  stress-straun  relation.  Using  the  velocity- 
stress  formulation  (Achenbach,  1975),  these  are  given  by 

d  d 

= /j  +  =  (1) 

y,,=  , . ]  12, 

. j,  ,3, 

where  Einstein’s  summation  convention  is  used.  J  is  the  dimension  of  the  problem, 
p  is  density,  and  A  and  p  are  Lame’s  parameters,  fj  are  body  forces  and  vj  and 
cTji  are  velocities  and  stresses  respectively. 

According  to  Tessmer  et.  al.  (1992)  and  Fornberg  (1988,  b).  a  2D  transfor¬ 
mation  of  a  rectangular  (^,  rj)-coordinate  system  into  a  curved  (x.  i)-coordinate 
system  is  a  computational  convenient  way  of  introducing  the  surface  topography. 
The  rect2mgular  (^,  7j)-system  is  limited  by  ^  =  0,  ^  ■  »?  =  0  and  rj  =  rjmaT 

(see  Fig.  1).  The  (x,  s)-sy8tem  is  curved  or  strticked  in  such  a  way  that  the  x-aods 
coincides  with  the  topography  at  the  top  of  the  system,  where  stretching  is  at  its 
maximum.  Below  the  surface  the  undulations  of  the  x-axis  decrease  linearly  with 
depth.  At  the  bottom  of  the  system,  the  x-axis  is  plane  (see  Fig.  2).  When  the 
positive  direction  along  the  vertical  axis  in  both  systems  is  upwards,  the  mapping 
from  the  (^ ,  »j)-system  to  the  (x,  j)-system  can  be  given  by 


II 

(4) 

(5) 

7ma» 
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where  »«  the  topography  function.  Applying  (4)  and  (5),  we  have,  for  a 
differentiable  function  f(x,z) 


dz  dij  dx ' 

U.  = 

dz  dr)  dz 


iFrom  expressions  (4)  and  (5)  we  get  (see  Appendix  1) 


BiO  = 


_  Vmax 
dz  -  zoiO' 


(6) 

(7) 


(8) 

(9) 

(10) 


Expanding  eqs.  (l)-(3)  by  the  chain  rule  and  setting  the  dimension  J  =  2 
leads  to  the  equations  in  the  medium  (see  Appendix  2) 

du  dCtx  .  -.dCxx  ,  /ii\ 

'ai  ”  di  dv  dv  *” 

=  ^  +  +  + 

^  =  (A  +  +  A5(0|^.  (13) 

d<Txz  ( dw  dw  du  ^ 

^=.(^  +  A«,,)^  +  S(0^),  (15) 

where  p  is  the  density  and  A  and  p  are  Lame’s  parameters,  fx  and  ft  are  the 
components  of  the  body  forces,  and  u  and  w  are  the  two  components  of  the  particle 
velocity.  <Txt,  <Zti  and  Cxt  are  the  stress  components. 
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Free  Surface  Boundary  Conditions 


The  free  boundzury  conditions  on  the  velocities  at  a  locally  horizontad  surface  read 


with  X  and  z  the  horizontal  and  vertical  coordinates  respectively. 

In  order  to  apply  these  conditions  at  a  topography  surface,  one  has  to  enforce 
them  at  eaich  point  on  the  surface  of  a  rotated,  local  coordinate  system  {x',z').  The 
orientation  of  the  system  (i',  z')  is  such  that  the  z'-axis  is  normal  to  the  surface 
at  each  point.  After  having  imposed  the  conditions  (16)  and  (17)  in  the  (x',z')- 
system,  they  have  to  be  rotated  back  to  the  (^,  7;)-system  by  the  transformation 

v  =  Av\  (18) 

where  v  and  v  '  are  the  particle  velocity  vectors  in  the  (^,Jj)-  and  the  {x',z')- 
systems  respectively.  A  is  the  rotation  matrix,  given  by 

A=(''“'‘  ‘“tV  (19) 

\-sin4>  cosp  J 

0  is  the  angle  of  rotation  between  the  actual  (r',  2')-system  and  the  (^,  77)-system. 
We  also  have  the  relation  tan0  =  52o(^)/5^.  Details  from  the  calculations  of  the 
rotation  aure  given  in  Appendix  3. 

We  arrive  at  the  conditions  for  the  free  surface  topography  given  in  the  com¬ 
putational  (^,  7/)-system  by 


Note  that  in  the  case  of  a  plane  surface,  i.  e.  zq(0  =  constant,  the  eqs.  (20)  and 
(21)  coincide  with  eqs.  (16)  and  (17). 

Numerical  Discretization 

Spatial  partial  differentiation  is  achieved  by  cost-optimized,  dispersion-bounded, 
high-order  finite  difference  operators  on  a  staggered  grid.  For  time  stepping  a  leap¬ 
frog  technique  is  used.  For  details  on  the  numerical  discretization  of  the  elastic 
equations,  the  reader  is  referred  to  Sguazzero  et.  al.  (1989).  To  avoid  artificial 
reflections  from  the  computational  boundaries,  the  multiplication  of  velocities  and 
stresses  by  exponentially  decreasing  terms  near  the  edges  works  satisfactorily.  For 
the  numerical  dispersion  relations,  the  stability  limit  and  bandwidth  introduced 
by  the  discretization,  the  reader  is  referred  to  Sguazzero  et.  al.  (1989). 

The  main  difference  in  the  discretization  of  eqs.  (20)  and  (21)  compared  to 
that  of  eqs.  (16)  and  (17)  is  tliat  we  now  stagger  the  vertical  velocity  component 
w  one  half  grid  length  downwards  instead  of  one  half  grid  length  upwards.  Modi¬ 
fications  according  to  this  change  have  to  be  taken  into  account  in  discretizing  the 
eqs.  (11)-(15)  inside  the  medium.  By  doing  this,  we  still  get  explicit  expressions 
for  the  variables  u,  w,  (Txx>  and  (Xxz  at  each  time  step  of  the  integration.  Sec¬ 
ond  order,  central,  uniform,  staggered  finite  difference  operators  (Fornberg,  1988, 
a)  are  retained  for  discretizing  the  free  boundary  conditions  eqs.  (20)  and  (21). 
Approaching  the  free  boundary  from  below,  we  apply  successively  lower  order  in 
the  central,  staggered  finite  difference  operators  used  for  discretizing  the  spatial 
derivatives.  Experiments  were  also  done  with  8th  order,  non-centred  finite  dif¬ 
ference  operators  according  to  Fornberg  (1988,  a)  for  discretization  near  the  free 
surface,  without  getting  stable  results  in  time.  A  pictorial  view  of  the  definition 
of  the  particle  velocities  on  the  numerical  grid  is  shown  in  Fig.  3.  The  staggered 
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deiiaitioQ  of  u  and  w,  the  horizontal  and  vertical  velocity  components  respectively, 
is  shown  here. 


Automatization  of  our  procedure  follows  from  the  fact  that  the  rotation  of 
each  local  coordinate  system  is  given  in  terms  of  derivatives.  Our  way  of  treating 
surface  topography  therefore  becomes  global.  One  does  not  need  to  consider  eewh 
spesific  surface  point  in  the  procedure.  As  also  the  medium  equations  are  solved 
numerically  by  the  same  procedure  for  the  whole  grid,  our  progreun  is  global  with 
regard  to  the  points  of  the  elastic  medium.  The  only  distinction  is  made  between 
the  surface  points  and  the  points  inside  the  medium. 

Another  superior  property  in  our  procedure  is  the  application  of  the  optimized 
finite  difference  operators  due  to  Sguazzero  et.  al.  (1989).  One  can  choose  a  certaun 
maximum  error  in  the  numerical  group  velocity  inside  the  medium,  where  the  op¬ 
erators  are  applied.  In  our  runs,  this  maximum  error  is  chosen  to  be  1.5  %.  Given 
a  chosen  error  for  the  numerical  group  velocity,  together  with  a  certain  chosen 
order  of  the  finite  difference  operator,  one  is  also  ensured  that  the  computational 
cost  of  the  resulting  operator  has  been  minimized. 


•  • 
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Numerical  Examples 

Results  of  numerical  examples  are  shown  in  Figs.  4-14.  The  two  particle  velocity 
components,  the  divergence  and  curl  are  displayed  in  snapshots  10  seconds  after 
a  pressure  wave  has  been  initiated  from  a  Gaussian  point  source  with  a  centred 
frequency  of  2.5  Hz.  The  source  is  placed  75  km  from  the  left  edge  of  the  grid 
and  2  km  below  the  surface.  The  total  size  of  the  computational  (^,  T;)-grid  is  150 
km  times  400  km,  with  751  grid  points  vertically  and  2001  points  horizontally. 
The  grid  length  in  each  direction  is  200  m.  Distances  of  the  snapshots  are  given 
in  kilometers.  In  each  snapshot,  only  the  upper  40  kra  of  a  100  km  horizontal 
subsection  is  shown  displaying  the  current  wave  front.  A  grey  tone  scale  is  given 
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for  each  variable.  Here,  divergences  and  curls  are  scaled  up  by  a  factor  of  10^ 
compaured  to  the  particle  velocities. 

The  surface  topography  used  in  our  runs  (Fig.  4)  is  generated  by  a  ID  real¬ 
ization  of  a  von  Karman  random  medium  (Charrette,  1991).  The  parameters  are 
order  =  1.0,  correlation  distance  =  10  km  and  standard  deviation  (RMS)  =  200 
m.  Distances  in  Fig.  4  are  given  in  kilometers. 

The  parameters  of  a  von  Karman  random  medium  can  be  understood  in  the 
following  way.  The  smoothest  version  of  the  medium  is  attained  for  order  =  1.0, 
while  the  strongest  variations  (self-similarity)  is  attained  for  order  =  0.0.  The 
correlation  distance  is  the  longest  distance  for  which  two  points  of  the  function 
in  question  are  highly  correlated.  Points  separated  by  more  than  the  correlation 
distance  are  close  to  uncorrelated.  The  standard  deviation  is  the  root  mean  square 
(RMS)  of  the  maximum  permitted  variation  from  the  average.  Also,  see  Frankel 
<k  Clayton  (1986). 

Too  abrupt  changes  in  the  topography  over  short  distances  leads  to  instability 
of  our  procedure.  Computing  the  derivatives  by  higher  order  operators  will  then 
lead  to  difficulties,  resulting  in  unwanted  oscillations  in  the  computations.  An 
experimental  measure  for  the  stability  limit  of  our  procedure  can  be  obtained  for 
each  type  of  topography  by  varying  parameters  characterizing  its  surface.  For 
our  von  Karman  curve,  experiments  were  done  by  varying  the  RMS,  while  keeping 
order  =  1.0  and  correlation  distance  =  10  km.  The  results  lead  to  an  approximate 
stability  limit  of  220  m  as  the  upper  bound  for  the  RMS.  As  mentioned  above, 
a  value  near  this  one  (200  m)  was  chosen  for  our  synthetics,  in  order  to  retain 
maximum  variations  in  the  topography. 

For  each  of  the  two  initial  runs,  a  different  type  of  medium  model  has  been 
applied.  Fig.  9  shows  the  first  model,  where  the  medium  P-velocities  are  defined 
by  different  linear  functions  vertically.  There  is  a  constant  low  velocity  layer  at 
the  upper  800  m,  where  the  P-velocity  is  4850  m/s.  In  the  next  800  m  below  this, 
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the  P-velocity  increases  linearly  from  4850  m/s  to  6150  m/s.  Then,  from  depth 
1.6  km  to  depth  35  km,  there  is  a  linear  increase  in  the  P-velocity  from  6150  m/s 
to  7000  m/s.  At  depth  35  km,  which  signifies  the  Moho,  we  have  a  discontinuity 
in  the  medium  P-velocity.  ^From  here  it  increases  linearly  from  8000  m/s  to  8500 
m/s,  which  is  the  P-velocity  at  the  bottom  of  the  grid,  150  km  below  the  top. 
The  snapshots  in  Figs.  5-8  have  been  generated  using  this  model. 

The  medium  model  applied  in  the  second  run  is  displayed  in  Fig.  14.  The  same 
linear  variations  in  the  medium  P-velocity  as  in  Fig.  9  are  applied,  but  additionally 
the  value  at  each  point  has  been  randomized  by  a  2D  von  Karman  realization 
(Charrette,  1991).  The  parameters  used  for  this  randomization  were  order  =  0.3, 
horizontal  correlation  distance  =  10  km,  vertical  correlation  distance  =  2.5  km 
and  RMS  =  3  %.  Thus,  an  apparent  anisotropy  has  been  included  in  the  medium 
by  letting  the  horizontal  correlation  distance  be  4  times  larger  than  the  vertical 
one.  The  snapshots  in  Figs.  10-13  have  been  generated  using  this  model.  Here  the 
waves  have  more  irregular  shapes.  In  each  of  the  runs  the  medium  shear  velocity 
is  defined  by  the  P-velocity  divided  by  \/3,  and  the  density  is  given  as  a  linear 
function  of  the  P-velocity. 

In  the  runs  of  this  paper  no  random  velocity  variations  were  applied  in  the 
upper  1600  m  of  the  medium.  Test  runs  with  random  velocity  variations  all  the 
way  up  to  the  free  surface  were  performed.  These  showed  scattered  surface  waves 
of  unrealistically  strong  amplitudes.  The  reason  for  this  may  partly  be  that  no  in¬ 
trinsic  attenuation  is  included  in  the  modeling.  Under  real  conditions  attenuation 
is  usually  quite  strong  near  the  surface  due  to  the  presence  of  a  weathering  layer. 

Snapshots  from  both  runs  reveal  mode  conversions,  back-scattering  and  asym¬ 
metric  P-to-P  and  P-to-S  reflections.  From  this  we  deduce  that  free  surface  to¬ 
pography  can  cause  these  effects.  Most  often  such  effects  are  discussed  only  in 
terms  of  inter-crust  heterogeneities.  Likewise,  the  fundamental  mode  Rayleigh 
(Rg)  wave  propagation  in  the  uppermost  part  of  the  crust  is  very  sensitive  to  the 
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roughness  of  the  local  topography,  as  synthesized  and  discussed  in  Ruud,  Husebye 
and  Hestholna  (1993). 

In  order  to  visualize  effects  of  including  surface  topography,  we  present  snap¬ 
shots  of  the  vertical  velocity  component  for  plane  surfaces  in  Figs.  15  and  16. 
Fig.  15  shows  this  with  the  medium  model  of  Fig.  9  applied,  and  Fig.  16  shows  it 
with  the  medium  model  of  Fig.  14  applied.  It  is  interesting  to  note  the  differences 
from  the  corresponding  snapshots  of  the  vertical  velocity  component  with  added 
topography  (Figs.  6  and  11).  In  Figs.  15  and  16  we  note  that  the  scattered  surface 
waves,  seen  as  a  black/white  banding  near  the  surface  in  Figs.  6  and  11,  are  almost 
absent,  particularly  in  the  horizontal  range  from  50  to  90  km. 

Fig.  17  display  seismograms  for  all  four  runs.  40  seconds  simulation  time 
is  used  to  include  adl  significant  wavelets  emanating  from  a  source  with  central 
frequency  of  2.5  Hz  at  depth  2  km.  In  each  case  6  receivers  are  located  at  the 
surface  starting  at  90  km  to  the  right  of  the  horizontal  position  of  the  source, 
with  a  spacing  of  800  m  between  them.  In  section  a)  of  Fig.  17  seismograms 
for  a  plane  surface  with  the  velocity  model  of  Fig.  9  are  shown  (no  topography 
or  medium  randomization).  Here  scattered  waves  are  absent,  and  the  following 
phases  are  clearly  seen:  Direct  P-wave  at  about  14  seconds,  P-wave  reflection 
from  the  Moho  at  17  seconds,  a  converted  P-to-S  reflection  from  the  Moho  at  23 
Ccconds  as  well  as  surface  waves  in  the  interval  26  to  37  seconds.  Here  the  wave 
with  highest  amplitude  and  lowest  frequency  is  the  fundamental  mode  Rg  wave, 
while  the  surface  waves  with  weaker  amplitudes  and  higher  frequencies  represent 
higher  order  modes  (Ruud  et.  al.,  1993). 

In  the  next  sections  of  Fig.  17  we  see  how  the  recordings  increasingly  get  dis¬ 
turbed  by  scattered  waves,  b)  shows  corresponding  seismograms  for  a  plane  surface 
with  the  velocity  model  of  Fig.  14  (no  topography  with  medium  randomization). 
Here  the  reflected  P-wave  from  the  Moho  has  an  abnormally  high  amplitude  on 
some  sensors.  This  is  probably  due  to  an  accidental  focussing  effect  of  the  velocity 


randomization.  The  two  last  sections  of  Fig.  17  (c)  aind  d))  contain  seismograms 
with  the  topography  of  Fig.  4  included,  c)  is  for  the  velocity  model  of  Fig.  9 
(topography  without  medium  randomization),  and  d)  is  for  the  velocity  model  of 
Fig.  14  (topography  amd  medium  randomization).  In  both  cases  the  fundaunental 
mode  Rg  wave  is  strongly  disturbed,  only  the  lowest  frequency  part  can  be  recog¬ 
nized.  Towards  the  end  of  the  two  sections  (35  to  40  seconds)  some  strong  coda 
waves  can  be  observed.  These  are  backscattered  surface  waves  moving  from  right 
to  left  in  the  model.  They  are  created  by  scattering  of  other  waves  at  the  free 
surface  topography.  In  d)  we  do  not  see  the  same  pattern  as  in  b)  from  the  acci¬ 
dental  focussing  effects  due  to  velocity  randomization.  This  is  because  the  random 
models  are  two  different  realizations  using  the  same  statistical  parameters. 

The  snapshots  shown  in  Figs.  5-8,  10-13,  15  and  16  were  arrived  at  after  4 
hours  and  20  minutes  on  the  IBM  3090  /  600  S  VF  computer*  located  at  IBM 
Bergen  Environmental  Sciences  and  Solutions  Centre,  Bergen,  Norway.  This  is 
the  location  for  adl  numerical  experiments  performed  for  this  work. 

Conclusions 

The  basic  assumption  of  the  calculations  is  that  the  topography  function  must  be 
smooth,  i.  e.  the  first  derivative  of  2o(0  must  exist.  If  the  function  varies  too  much 
over  a  certain  region,  Gibb’s  phenomenon  will  arise,  and  the  resulting  oscillations 
will  destroy  the  simulation.  The  limit  for  this  to  occur  is  a  matter  of  experimen¬ 
tation.  Given  that  this  stability  condition  is  fulfilled,  the  above  procedure  is  a 
powerful  tool  for  the  finite  difference  modeling  of  an  elastic  medium  including  free 
surface  topography. 


•  • 
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Appendix  1 


For  the  equations  in  the  medium  we  need  |^ .  and  .  They  aure  found  from 


dx  dxdT)_ 

dx  ^  drfdx 
dz  —  0 

dx  ^  dr)  dx 


dx  d^  _  f) 

d^  dz  ^  dr)  dz 
dz  di  ^^7  _  j 

dlTz'^^Tz  ~ 


This  leads  to 


dx 

dx 


--/J 
d^'  ' 


dz 

dz 


with 


J  = 


dx  dz  dx  dz 


d^  dr)  dr)  d^ 

With  our  choice  of  mapping  functions,  eqs.  (4)  and  (5),  we  get 


dx 

=  1. 

dx 

dr) 

dz 
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Ci3 

U 
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rn 

dz 

r)max 

dr) 

0. 

zo(0 

r)max 
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J  _  ^  _  ^o(^) 

dr)  r)mas 

^From  this  we  get  the  expressions  eqs.  (8)-(10). 
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Appendix  2 


Applying  the  chain  rule  on  eqs.  (I)-(3)  with  the  dimension  J  =  2  leauls  to 


^  di]  dr) 

^  dt  dx  dr)  dx  d^  dz  dt)  dz  *  ’ 

dw  _  dvx,  d^  d<Tzt  dr)  diXttd^  d<T,t  dr) 

^  dt  d^  dx^  dr)  dx^  d^  dz^  dr)  dz^ 


dt 

d<T,z 

dt 

dOxz 

dt 


du  d^  du  dr)\  / dw  d^  dw  dr) 

d^  dx  ^  dr)  dx )  ^  V  dz  dr)  dz 


=  A 

=  /i 


( dud^  dudr)\ 

/  dw  d^  dw  dr)  du  d^ 

\  di  dr)  dx  ^  d^  dz^ 


( . 

\d^dz^ 

du  dr)\ 
dr)dz) 


dw  dr) 
dr)  dz 


)■ 

)■ 


Substituting  for  d^/dx,  d^/dz,  dr)/dx  and  dr)ldz  from  eqs.  (8)-(10),  we  get  the 
medium  equations  eqs.  (11)-(15). 
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Appendix  3 

Suppose  a  vector  v  is  given  in  a  coordinate  system  with  basis  vectors  i  and  j.  This 
system  is  then  rotated  an  angle  0  into  a  new  coordinate  system,  with  basis  vectors 
i  '  and  j  In  this  new  system  the  vector  v  is  denoted  by  t?'.  We  then  have  the 
connections 

i  '  =  cos  •  i  +  sin  ^  •  j 

j  '  =  —  sin  i  +  cos  ^  •  j 

with  I  being  the  vector  Correspondingly,  we  have 

=  CO80  •  i  '  -  sin<^  •  J  '  1  _  .  _i  r, 

=  sin  0  •  J '  +  cos  <t>  ■  j  }  ~  ’ 

with  I '  being  the  vector  (i  ',  j  ')  amd  A~*  being  the  inverse  of  the  rotation  matrix. 
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The  coordinate  transformation  for  v  is  described  by  v  =  Av  ' .  Hence,  v  '  = 
or,  with  u'  and  w'  the  components  of  v and  u  and  w  the  components  of 


u'  =  cos  ^  •  u  —  sin  ^  •  u;, 
u;'  =  sin  ^  •  u  +  cos  ^  •  w. 


Applying  the  chain  rule  on  a  differentiable  function  f,  we  have 


dx'  34  dx'  dr)  dx‘ 

=  —  cos.^+  -^sin^, 

3^  dt) 

3z'  3^  3z‘  3v  3z‘ 

f  f 

=  ^(-sin^)  +  ^cos<^. 
5^  3t} 


We  have  to  enforce  the  free  boundary  conditions  on  the  velocities  into  the  (x',z')- 
system,  where  the  r'-axis  is  normal  to  the  surface  at  the  local  point,  i.  e. 


A  3u' 
X  +  2fi  3x' 


If  the  chain  rule  is  applied  as  above,  we  get 


3u'  .  du'  3w'  3w'  . 

•j^(- sm «  +  CO, ♦  -  «u 


3w'  .  3w'  X 
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/  3u'  3u'  .  A 
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Now  we  use  the  expressions  given  above  for  ti'  and  u;'  which  we  got  from  the 
rotation.  Then  we  obtain 

•  2  .  2  , 

-8in0co8  0—  +sm'‘0-^  +  co6‘‘0-T - sm^cos^— 

a<  Of]  at) 

■  ,  ,du  ^  ,dw  ■  i ,du  .  ,  dw 

=  -sm^cos^—  -  cos'*  —  sm  <i>- - sm(>cos0-r-. 

dr)  at) 


.  .du  .  ^  dw  du  ■,  dw 

-sm  <p-^  ~  sm<t>cos<p-^  +  sm  ^cos^—  +  cos' 

a<  d^  dr)  dr) 


f  a  ■  a  .  •  a 

<  cos  0-^^  -  sm 0 cos ^ +  sm  4> cc»4>—  —  sm*  (j> 


X  +  2fi 


dr) 


dr)] 


The  first  of  these  equations  immediately  leads  to 


du  _  dw 
dr)  ~  d^ 


The  second  equation  leads  to 
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dw 

dr) 


=  tan<^(  1  + 


(■ 


X  +  2fiJ  d^ 


tan  ^  (  1  + 


A  \du 


X  +  2)1  J  dr)’ 


•  •  • 


•  • 


after  dividing  by  cos*  <t)  and  restructuring.  This  division  requires  the  assump¬ 
tion  (p  #  d:ir/2,  which  means  that  the  topography  cannot  have  vertical  subsec¬ 
tions,  i.  e.  the  topography  function  must  be  single-valued.  Given  the  condition  of 
smoothness  for  zo(?)i  this  is  a  reasonable  assumption.  Now  we  apply  the  first  of 
these  equations  into  the  second  one  and  use  the  equality  tan^  =  dzo{^)/d^.  Then 
we  arrive  at  eqs.  (20)  and  (21)  as  a  possible  set  of  free  boundary  conditions  at  a 
topography  surface. 
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1th  velocity  component  u  defined  as  *. 
2nd  velocity  component  w  defined  as  =. 


Figure  3;  Particle  velocity  definition  on  the  numerical  grid. 
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Figure  4:  Surface  topography  used  in  the  simulations. 


Horizontal  velocity 


Figure  5:  Snapshot  of  horizontal  particle  velocity  10  seconds  after  a  pressure 
wave  is  initiated  from  a  Gaussian  point  source  with  depth  2  km  below  the 
surface  topography.  The  medium  model  is  displayed  in  Fig.  9. 
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Figure  11:  Snapshot  of  vertical  particle  velocity  for  the  same  situation  as 
in  Fig.  10. 
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Fisnrp  14:  Random  medium  model  of  pressure  velocity  used  for  the  stiap- 
siiots  in  Figs.  10  -  Id. 
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Figure  17:  Seismograms  for  all  four  runs.  In  each  section  six  traces  with 
800  m  separation  are  given,  starting  at  90  km  to  the  right  of  the  horizontal 
position  of  the  source.  They  display  the  vertical  particle  velocity  at  the 
free  surface.  The  four  sections  represent:  a)  plane  surface  with  no  random 
velocity  realization,  b)  plane  surface  with  random  velocity  realization,  c) 
topography  surface  with  no  random  velocity  realization  and  d)  topography 
surface  with  random  velocity  realization. 
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